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THEORY  AND  NUMERICAL  MODELING  OF 
LOW-FREQUENCY  ACOUSTIC 

SCATTERING  FROM  BUBBLE  PLUMES  NEAR  THE  SEA  SURFACE 

1.  INTRODUCTION 

In  this  report,  we  develop  the  physical  theory  and  numerical  techniques  for  modeling  low- 
frequency  scattering  from  acoustically  penetrable  bubble  plumes,  including  surface-image  and  res¬ 
onance  effects.  The  model  plumes  can  have  arbitrary  shapes,  interior  sound  speed  fields,  and 
depths.  In  addition,  a  pressure  relief  boundary  condition  can  be  invoked  at  the  air/sea  interface, 
which  is  taken  to  be  flat. 

1.1  Background 

At  low  sea  states  and  frequencies,  acoustic  sea-surface  reverberation  is  fundamentally  just 
rough-surface  scattering.  However,  there  is  now  ample  evidence  that  with  increasing  sea  state 
and/or  frequency,  a  second  scattering  process  appears  and  eventually  dominates.  Reference  1 
provides  a  complete  review  of  the  subject  and  we  will  not  go  further  into  the  matter  here  except  to 
say  that  it  is  now  virtually  certain  that  subsurface  bubble  structures  are  the  scattering  mechanism 
responsible.  Even  with  that  established,  the  acoustic  situation  remains  somewhat  unclear  because 
scattering  from  such  structures,  particularly  if  they  are  resonant,  is  not  well  understood. 

In  the  present  context,  “low  frequency”  simply  means  well  below  the  resonant  frequencies  of 
the  individual  bubbles  (e.g.,  /  <  900  Hz,  for  10-micron  bubbles).  In  this  regime,  small  concentra¬ 
tions  of  bubbles  in  the  water  (e.g.,  less  than  1  %  air  by  volume)  have  a  negligible  impact  on  the 
density  but  produce  a  dramatic  increase  in  the  medium’s  bulk  compressibility  (a  sensitive  function 
of  the  air  volume  fraction)  and  a  commensurate  reduction  in  sound  speed  [2,  3]. 

The  work  presented  here  was  undertaken  to  clarify  the  scattering  properties  of  one  type 
of  subsurface  bubble  structure,  the  so-called  “intermediate  bubble  plume.”  These  are  localized 
collections  of  air  bubbles  that  get  injected  into  the  surface  scattering  zone  by  plunging  or  breaking 
waves.  They  are  characterized  by  void  fractions  on  the  order  of  1%  [4,  5]  so  that  the  interior 
sound  speed  may  be  as  much  as  five  times  smaller  than  the  sound  speed  outside.  Because  of 
these  remarkably  low  interior  sound  speeds,  intermediate  bubble  plumes  can  function  as  resonant 
cavities  even  at  low  frequencies.  Only  a  few  short-lived  clouds  of  this  type  might  be  needed  to 
produce  a  significant  contribution  to  the  scattering  cross-section  from  near  the  ocean  surface. 

1.2  Approach 

The  essentials  of  our  approach  are  all  present  in  the  standard  treatment  of  oiie-dimensional 
scattering  from  a  uniform  object;  i.e.,  a  square  well.  Appendix  A  reviews  the  familiar  solution 
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method:  (a)  write  the  fields  inside  and  outside  the  scatterer  as  suitably  general  parameterized 
expressions,  (b)  impose  physical  continuity  at  the  boundaries,  (c)  solve  the  boundary  continuity 
conditions  for  the  unknown  parameters,  and  (d)  insert  these  parameter  values  into  the  expression 
for  the  outgoing  field  far  from  the  scatterer.  Figure  1  shows  the  result  of  such  a  calculation  for  a 
well  with  a  width  of  3  m,  exterior  sound  speed  of  1500  m/s  and  interior  sound  speed  of  100  m/s. 
The  transmission  coefficient  is  a  smooth  function  of  frequency  punctuated  by  peaks  where  the 
internal  resonances  occur.  The  narrowness  of  these  resonance  peaks  (i.e.,  the  resonator’s  Q-factor) 
depends  on  the  size  of  the  object  (i.e.,  the  width  of  the  well)  and  the  sound  speeds  involved. 


Fig.  1  -  Magnitude  of  the  trauismission  coefficient  for  a  uniform  1-D  object.  Sound  speeds: 
exterior  1500  m/s,  interior  100  m/s;  object  size,  3  m.  Dashed  vertic^ll  lines  Me  spMed  at 
A/  =  0.5  X  (interior  speed) /(object  size).  See  Appendix  A  for  details. 


The  scattering  problem  in  this  report  is  much  more  complicated  than  that,  but  this  complex¬ 
ity  is  really  not  fundamental:  it  is  merely  a  byproduct  of  dealing  in  higher  dimensions  without  the 
simplifying  symmetry.  In  fact,  the  solution  follows  the  same  basic  pattern.  General  expressions 
for  the  interior  and  exterior  fields  are  given  in  terms  of  the  Green’s  functions  for  the  two  regions. 
Physical  continuity  at  the  boundary  is  embodied  in  a  pair  of  integral  equations  over  the  plume’s 
surface.  The  solution  of  these  boundary  integral  equations  yields  distributions  of  surface  sources 
that,  in  effect,  reradiate  the  scattered  field.  Finally,  the  far-field  scattering  amplitude  appears  as  a 
surface  integral  over  these  equivalent  distributions.  The  resulting  spectra  are  qualitatively  similar 
to  Fig.  1,  although  the  resonances  are  not  nearly  so  sharp,  regularly  spaced,  or  uniform  in  height 
for  realistic  bubble  plumes. 

Our  statement  of  the  problem  for  this  report  has  its  limitations.  Among  these  are  the  neglect 
of  curvature  and  roughness  in  the  sea  surface  (it  is  taken  to  be  flat  and  smooth),  time  evolution 
of  the  plume  itself,  and  very  large  interior  sound  speed  gradients.  By  omitting  these  effects  from 
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the  current  formulation  we  do  not  mean  to  dismiss  them  as  unimportant.  We  are  simply  bounding 
the  problem  by  deferring  these  components  for  consideration  at  a  later  time. 

2.  INTEGRAL  FORMULATION 

The  plume’s  boundary  is  a  surface  <S  with  outward  normal  direction  n.  The  exterior  region 
£  has  a  constant  sound  speed  and  density.  The  interior  I  also  has  a  constant  density  but  its  sound 
speed  is  a  function  of  the  air  volume  fraction  and  thus  may  vary  from  point  to  point.  Since  the 
sound  speed  throughout  the  bubbly  water  in  X  is  lower  than  that  of  the  normal  water  in  £  (Section 
6),  the  opposite  is  true  of  the  wavenumber  k  =  uj/c:  its  smallest  values  are  found  in  the  exterior. 
The  situation  is  sketched  in  Figs.  2  and  3. 


Fig.  2  -  Sketch  of  a  bubble  plume  showing  its  interior  Z,  exterior  £,  surface  S,  outward 
surface  normal  n(r),  and  linear  dimension  D 


If  an  actual  sound  speed  discontinuity  exists  between  bubbly  and  normal  water,  as  shown  on 
the  right-hand  side  of  Fig.  3  (a  and  b),  that  discontinuity  marks  the  surface  of  the  plume.  If  not. 
S  may  be  taken  to  lie  anywhere  beyond  the  limit  where  c  is  effectively  constant.  (For  that  reason, 
the  dotted  line  on  the  left-hand  side  of  the  figure  could  just  as  well  be  drawn  farther  to  the  left.) 
This  report  treats  both  the  “sharp”  and  “fuzzy”  boundary  cases. 

2.1  Equivalent  Source  Distributions 

Kittappa  and  Kleinman  [6]  have  used  classical  methods  to  formulate  the  boundary  integral 
equations  that  govern  continuous  wave  (CW)  scattering  from  a  penetrable  body  of  very  general 
shape.  They  allow  the  surface  to  have  creases  and  corners,  but  restrict  their  attention  to  uniform 
homogeneous  interiors.  We,  however,  will  focus  on  the  contrasting  case  of  a  smooth  surface  and  an 
inhomogeneous  interior,  since  this  more  closely  approximates  the  reality  of  bubble  plumes  in  the 
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Fig.  3  -  Sketches  of  sound  speed  and  wavenumber  as  functions  of  one  of  the  spatial  coordi¬ 
nates.  The  dotted  lines  indicate  the  boundary  surface  S. 


ocean.  By  a  straightforward  adaptation  of  their  approach  to  this  case,  one  arrives  at  two  integral 
equations,  valid  over  the  surface  of  the  scatterer,  that  relate  the  total  acoustic  velocity  potential 
U  and  its  normal  derivative  V  to  their  corresponding  incident  values  Vine,  Vine-  These  boundary 


integral  equations  are  [6,  Eqs.  (3.5)  and  (3.19)] 

f/inc(y)  + j^d5(x)(A(x,y)17(x)-B(x,y)F(x))  =  {p/pint)U{y)  (la) 

and 

ViTic(y)  + jf  d5(x)(C(x,y)(7(x)-D(x,y)F(x))  =  {ptPext)V{y)  (lb) 

for  y  E  S,  where  U{y)  and  V(y)  represent  limiting  values  as  y  — ►  5  from  within  £,  the  average  of 
the  two  densities  is  denoted  p  —  {pint  +  Peit)/2,  and  the  kernels  are 

^(x,y)  =  ^^^{Gext(x,y)  -  (Pext/Pmt)  Gj„f(x,y)}  (2a) 

B(x,y)  =  {Geit(x,y)  -  G’i„t(x,y)}  (2b) 

I>{x,y)  =  -Q;^^{Gext{x,y)  -  {Pint/Pext)  Gint[x,y)}  (2d) 

with  ‘9/5n’  indicating  a  derivative  in  the  outward  normal  direction. 


As  observed  in  Ref.  6,  the  B  kernel  is  always  regular  while  C  is  weakly  singular,  i.e.,  it 
diverges  as  e  =  |x  -  y|  — >  0  but  in  a  milder  way  than  and  A  and  D  are  at  worst  weakly 
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singular  too.  (In  the  familiar  case  where  pint  =  Pext,  ^  and  D  are  regular.)  Since  the  differential 
area  element  dS  is  proportional  to  the  integrals  of  C,  A,  and  D  are  analytically  well-behaved, 
though  they  are  challenging  numerically. 

The  field  in  the  exterior  region  is  the  sum  of  incident  and  scattered  parts,  U  =  Ut„c  +  Usca, 
the  scattered  part  being  [6,  Eq.  (3.3)  evaluated  in  the  exterior] 

Uscair)  =  dSix)  idGext{^,r)/dnix)  U (x)  -  Gext{x.r)V (x))  (3) 

for  r  €  5,  in  terms  of  the  U,  V  surface  data  from  Eq.  (1).  In  this  expression,  the  surface  distribu¬ 
tions  U  and  V  appear  as  monopole  and  dipole  sources  for  the  scattered  field  and  are  “equivalent 
source  distributions”  in  that  sense. 

The  field  C/i„c  has  thus  far  been  left  unspecified.  Hereafter  we  focus  our  attention  on  the 
incident  field  from  a  point  source  at  s.  With  this  choice  and  with  the  s  dependence  made  explicit 
by  writing  f7(y|s)  and  y(y|s),  Eqs.  (1)  and  (3)  amount  to  the  matrix  equations 


{p/Pint)U{y\s) 
{p/PextW  {y\s) 


(  d^ffxl  (  S(x,y)  \  /  I/(x|s)  \ 

Is  (  -C(x,y)  D(x,y)  )  (  V{x|s)  )  "  ( 


for  surface  points  y  €  5  and 


Uscair\s)  =  l^dS{x)(  dGextix,r)/dnix),  -Gext(x,r)  ^  (5) 

for  exterior  points  r  €  S.  For  any  exterior  source  point  s,  these  equations  are  solvable  for  f7sca(r|s) 
at  any  exterior  field  point  r.  In  the  next  section,  we  concentrate  on  solving  them  in  the  far-field. 


2.2  Far-Field  Limit 

We  now  pass  to  the  far-field  limit.  With  the  origin  in  the  vicinity  of  the  scatterer,  we  allow  s 
and  r  to  recede  to  distant  reaches  of  £  so  that  s,  r  »  D.  Then  we  can  expand  the  Green’s  function 


Gext(y,s)  = 


exp{ffcextiy  -  s|} 
47r|y  —  s| 


and  its  normal  derivative  at  a  surface  point  y  to  first  order  in  (D/s)  to  obtain 


Gextiy,s) 

Gextiy,s) 


exp{ikexts} 

4ns 

exp{ikexts} 


Pis,y) 

Q(s,y) 


where 


-P(s.y)  =  exp{-tA:sx<sy} 


Q(s,y)  =  n(y)-s  P(s,y). 
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Physically,  this  just  means  that  the  spherical  incident  field  is  approximately  a  plane  wave  through¬ 
out  the  scattering  region.  The  phase  of  this  plane  wave  is  kext-y,  where  the  source  direction  is 
s  =  —kext-  This  estimate  is  valid  only  in  the  source  terms  on  the  right-hand  side  of  Eq.  (4). 
Similarly,  for  a  distant  receiver  and  a  surface  point  x,  the  estimates 

Geit(x,r)  ss  -^^^^^P(f,x)  (9a) 

and 

S^G,xt(x,r)  «  - — - Q(r,x)  (9b) 

apply  in  the  scattered  field  integral,  Eq.  (5).  As  a  result,  in  the  far-field  regime,  the  scattered  field 
has  the  asymptotic  form  [7,  8] 

.  (10) 

This  is  an  outgoing  spherical  wave  modulated  by  a  complex  amplitude  F(f  |kcj;()  —  the  scattering 
amplitude  or  “far-field  radiation  pattern”  in  the  direction  f  that  is  produced  by  incident  plane 
wave  radiation  with  unit  amplitude  and  wavevector  kext-  Figure  4  shows  that  k^xt  =  fceitkcit  = 
(w/cext)(— s).  It  follows  that  F(f  |kext)  is  the  amplitude  for  scattering  at  frequency  lj  from  direction 
— s  into  direction  f .  Because  the  frequency  is  fixed,  we  are  free  to  leave  the  kext  dependence  implicit 
and  write  the  scattering  amplitude  as  simply  F(f|  -  s).  This  single-frequency  far-field  scattering 
amplitude  is  what  we  will  model. 


Fig.  4  -  Sketch  of  the  quantities  involved  in  the  scattering  function 
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In  a  similar  way,  t/(y|8)  and  V(y|8)  have  the  far-field  asymptotic  forms 


and 


V"(yis) 


exp{27i:eit^} 

47rs 


(lla) 

(llb) 


where  t/(y|s)  and  V’(y|8)  are  factors  that  depend  only  on  the  direction  of  s.  In  the  far-field  limit, 
neglecting  spherical  spreading  factors  exp{ifceitr}/47rr  and  exp{ikexts}/^ns,  our  task  is  first  to 
solve  the  boundary  continuity  condition 


/  {p/pint)U {y\s)  \  f  /  -^(x,y)  B{x,y)  \(  17(x|s) 

\  (p/pext)V(yls)  j  Js  Y  -C:'(x,y)  £)(x,y)  j\^K(x|s) 


P{s,y)  \ 
Q(s-y) ) 


(12) 


obtaining  l/(x|s)  and  V(x|s)  throughout  x  €  «S,  and  then  to  use  these  in  calculating  the  scattering 
amplitude  according  to 

=  ii Is  (  nS)  )  ■ 

As  usual  in  scattering  work,  we  are  treating  the  exterior  region  as  if  the  local  value  of  c^xt 
prevailing  near  the  scatterer  actually  extended  to  infinity  in  all  directions.  In  environments  w'here 
that  approximation  is  valid,  the  scattered  field  is  expressed  simply  and  directly  in  terms  of  this 
scattering  amplitude.  Even  when  the  exterior  region  is  considerably  more  complex  (in  a  shallow 
water  waveguide,  for  example)  the  plane  wave  scattering  amplitude  remains  a  useful  theoretical 
and  computational  intermediate  quantity  [9,  10].  The  remainder  of  this  section  introduces  the 
modification  required  to  accommodate  the  simple  departure  from  an  infinite  homogeneous  exterior 
that  will  be  needed  for  scattering  from  a  near-surface  plume  —  namely,  the  presence  of  the  air/sea 
boundary. 


2.3  Air-Sea  Boundary  Condition 


We  will  take  the  flat  air /sea  interface  z  =  0  to  be  a  pressure- release  boundary.  Temporarily 
shifting  the  origin  up  from  the  center  of  the  plume  to  the  interface,  we  resolve  vectors  into  com¬ 
ponents  parallel  and  perpendicular  to  that  boundary,  e.  g.,  x  =  xy  -I-  xj_,  and  denote  the  image 
under  reflection  by  x  —  xy  -  xj^.  The  pressure  relief  condition  can  be  implemented  at  this  stage 
of  the  problem  simply  by  using  the  appropriate  half-space  exterior  Green’s  function  instead  of 
Eq.  (6).  All  of  the  above  results  remain  valid  if  we  simply  replace  Geit(x,y),  PCr.x),  and  Q(r,x) 
by  Gext{^,y)  —  Geit(x,y),  P(f,x)  -  P(f,x),  and  Q(f,  x)  -  Q(f,  x).  Explicitly,  this  means  using 

^  exp(ifcext|x-y|)  exp  (ifcext|x  -  y|) 

LreitlX,y)  =  - - j - j - — j - - -  ,  (14) 

47r|x-y|  47rlx-y| 

P(r,x)  =  -2isin(A:eitx  fx)exp{-ifcej,tX-fy}  ,  (15a) 

and 

Q(r.x)  =  -“^kext  (sin(A:ea.tX-fx)fy  -|- i  COS  (fceit  X-fx)  Tx)  •n(x)  exp{-?fre3-t  X-Fy  }  .  (15b) 

This  is  sufficient  if  the  interface  lies  entirely  in  the  exterior,  i.e.,  if  the  plume  does  not  break  the 
surface.  We  will  assume  so  and  proceed  to  treat  plumes  of  that  type. 
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3.  GREEN’S  FUNCTIONS 

In  this  section,  we  develop  expressions  for  the  Green’s  functions  and  their  gradients  that  are 
needed  in  the  A,  B,  C,  and  D  kernels  of  Eq.  (2).  For  this,  we  define  the  relative  vector  u  =  x  —  y 
and  exploit  the  fundamental  gradients  of  its  length  u\  thus 

u  , 

=  -u  , 


=  —^/h  ,  (ific) 

where  u  =  u/u  is  ttie  direction  of  u,  11  =  uu  is  the  projection  operator  in  that  direction*,  and 
^  =  1  —  n  is  the  complementary  projector  onto  the  plane  orthogonal  to  u. 

3.1  Exterior 


(16a) 

(16b) 


and 


a 

m 

a 

Sy 


u 


u 


a 


a 


Differentiation  of  Eq.  (14)  produces  the  following  expressions  for  the  half-space  exterior 
Green’s  function  and  its  spatial  gradients: 


and 


^  m  G'eit(x,y)  - 


G.xt(x,y)  =  ^  , 

4tt  \w  V  ) 

SEG.«(x.y)-  “'ir'  “V  ’ 

gy  Gext(x,y)  -  u-^e  UtJ  , 

_  ( p(tu)uu  -  q{w)l  p(u)u^U|  -  q{v)i 


47r 


g.W  _ 


(17a) 

(17b) 

(17c) 

(17d) 


in  ternr.3  of  the  phase  factors 


w  =  Areitlx  -  y| 

V  =  fceitjx  -  yl  , 


the  complex  functions 


and  the  unit  vectors 


g(x)  =  *x  - 1 
P(x)  =  X^  +  3g(x)  , 


Ut 


x-y 

Ix-yl 

x-y 

|x-y| 


(whose  vertical  components  point  downward  and  upward  as  indicated  by  the  arrow  subscripts)  as 
well  as  the  identity  and  reflection  operators  1  and  1  =  1  —  20^62  for  which  lx  =  x  and  ix  =  x. 
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3.2  Interior 


Assuming  small  interior  sound  speed  gradients  in  the  sense  that  |AA/A|  1;  i.e., 


ewKB  =  kintC**) 


,dkint{r) 


1  , 


we  may  approximate  the  interior  Green’s  function  by  the  WKB-type  form 


G'ine(x,y)  = 


e«V(x,y) 


The  required  gradients  are  then 


(19d) 


in  terms  of  the  spatial  derivatives  of  the  phase  term  </?.  With  the  line  segment  from  y  to  x 
parameterized  by 

z(0  =  y  +  ^u  •••  o<^<u,  (20) 

this  phase  can  be  expressed  as  the  line  integral 

<^(x^y)=/  d^kintiziO)  (21) 

Jo 

from  z(0)  =  y  to  z(u)  =  x.  Its  first-order  gradients  are  consequently 

—  =  -I-  fci„t(x)u  (22a) 

dip  dip 

i  =  •‘-I 

and,  with  second-order  spatial  derivatives  of  the  wavenumber  neglected,  the  second-order  gradient 
is 

as  -  ( u{f-^}  -h  (f -f  (f-u  -  fc,„t(x))^ )  -  ug(x)  -I-  2(g(x)-u)n  (22c) 

cfycfx  u 

in  terms  of  the  interior  wavenumber  gradient 


g(r)  = 


dkintir) 


and  its  integrals 


e  =  /  d^g(z(0) 

Jo 

f=-  r  d^^gizio)- 

u  Jo 
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Thus,  the  phase  and  its  gradients  can  be  computed  by  numerical  quadrature  and  then  used  to 
compute  the  interior  Green’s  function  and  its  gradients. 

With  these  explicit  forms  for  the  two  Green’s  functions,  we  now  have  all  of  the  components 
essential  to  the  analytic  specification  of  the  scattering  problem;  namely,  the  kernels  A,  B,C,  and 
D  in  the  continuity  condition  (Eq.  (12))  and  the  P  and  Q  functions  that  appear  in  Eq.  (12)  as 
inhomogeneities  and  also  as  quadrature  weights  in  the  scattering  amplitude  equation,  Eq.  (13). 
The  remaining  task  is  to  devise  a  means  of  getting  a  numerical  solution. 

4.  WIREFRAME  SOLUTION 


4.1  Continuity  Condition 


To  solve  the  continuity  condition  (Eq.  (12))  numerically,  we  first  ‘tile’  the  plume  surface.  This 
means  fitting  S  with  a  “wireframe”  of  N  plane  facets  (the  tiles)  that  are  approximately  uniform 
equilateral  triangles  (Fig.  5,  for  example)  with  centers  at  the  points  Xn.  We  then  approximate  the 
surface  integral  as  a  sum  over  the  wireframe,  thus  converting  Eq.  (12)  into  the  algebraic  expression 


U(s)  \  /  P(s)  \ 

V(s)  )  [  Q(s)  )  ■ 

Here  I  is  the  iV  x  AT  identity  matrix;  U(s),V(s)  and  P(s),Q(s)  are  the  N  x  I  column  matrices 


{p/pint)^  A  B 

(p/Pext)l  +  D 


(25) 


and 


and  A,  B,  C,  and  D  are  iV  x  iV  matrices  with  off-diagonal  elements 

Afnn  —  ■A(x,i,  Xm)  Sm  , 

Bmn 


Un{s) 

=  t;(Xn  |s)  , 

(26a) 

Vn{s) 

II 

(26b) 

Pn{s) 

=  P(s|x„)  , 

(26c) 

Qn{s) 

=  Q(slx„)  ; 

(26d) 

—  ■5(x^,  Xm)  Sjji  , 

=  C’(x^,  Xm)  *Sm  , 


(27a) 

(27b) 

(27c) 


and 


Dmn  —  U(Xn,  Xtji)  Sm, 


(27d) 


where  Sm  is  the  area  of  the  patch  of  the  plume’s  surface  subtended  by  the  mth  tile  (see  Appendix 
B).  These  m  ^  n  matrix  elements  can  be  calculated  directly  from  the  definitions  of  Exj.  (2) 
using  the  weas  and  normals  of  the  surface  patches  and  the  Green’s  function  gradients  that  were 
develop>ed  in  the  preceding  section.  The  m  =  n  elements  require  separate  consideration. 

4.2  Diagonal  Matrix  Elements 

The  diagonal  elements  are  complicated  by  the  weak  singularity  of  some  of  the  kernels  in 
Eq.  (12),  i.e.,  by  their  divergence  as  x  — ►  y  [8].  Expressions  for  them  are  derived  in  Appendix 
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C  and  axe  summarized  below.  In  these,  11^  and  Zm  are  the  mth  surface  patch’s  normal  vector 
and  center  depth,  kint,m  =  fctnt(xm)  is  the  interior  wavenumber  at  the  patch  center,  Wm  is  the 
integral  /  d5(x)l/|x  —  Xm|  over  the  mth  patch  (Appendix  B)  and,  as  before,  —  1. 

Also  appearing:  pm  =  “i-ZmKxti  which  measures  the  depth  of  the  patch  in  units  of  the  exterior 
wavelength;  =  e2-n„j,  which  represents  the  vertical  tilt  of  the  normal;  and  Pm  =  “m  gCxm), 
the  normal  component  of  the  wavenumber  gradient. 

The  diagonal  matrix  elements  are 


^mm  — 

(t) 

(28a) 

Bmm  — 

-  {•(feni.m  -  k„t)  +  (-^)  e"^}  (^) 

(28b) 

Cmm  — 

-  e.)  (^)  -  (^) 

(28c) 

Dmm  — 

(28d) 

The  terms  containing  embody  the  surface-image  effects.  With  the  matrix  elements  provided 
by  Eqs.  (26)  through  (28),  Eq.  (25)  can  be  solved  for  the  values  of  the  source  densities  U,V  on 
each  of  the  surface  tiles. 

4.3  Scattering  Amplitude  Equation 

After  Eq.  (25)  is  solved  for  U(s),  V(s),  all  that  remains  is  to  apply  Eq.  (13)  in  the  wireframe 
approximation;  i.e., 
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F(f|-s)  =  j;5„(Q(f,x„),-P(f,x„))  ^  ^ 

=  Q^(f)U(s)  -  P'^(f)V(s)  (29) 

where  Q  and  P  correspond  to  Q  and  P  with  components  weighted  by  the  surface  area  elements 
(Qn  =  SnQn,  Pn  =  S„Pn)  and  ‘  ^  ’  indicates  a  matrix  transpose. 

5.  MODELING  METHODS 

This  section  describes  the  numerical  methods  that  we  have  used  to  implement  the  scattering 
formulation  described  above.  We  should  point  out  now  that  the  pressure-relief  boundary  condition 
that  we  have  taken  such  pains  to  impose  at  the  sea  surface  may  be  omitted,  if  desired,  by  simply 
omitting  the  v,fim  terms  in  Eqs.  (17)  and  (28)  and  reverting  from  Eq.  (15)  back  to  the  elementary 
forms  in  Eq.  (8).  Indeed,  during  prograim  development  we  supported  this  choice  as  a  diagnostic 
option.  We  have  retained  that  option  so  that  the  final  model,  known  as  BIRPS  (Boundary  Integral 
Resonant  Plume  Scatter),  can  also  be  used  to  model  “free-space”  scattering. 

5.1  Construction  of  the  Wireframe 

All  the  wireframes  are  initially  spherical.  Their  construction  is  handled  by  an  algorithm  that 
begins  with  2"*  tiles  (m  >  2)  around  the  sphere’s  equator  and  repeatedly  halves  this  number  at 
higher  latitudes  to  control  the  tile  size,  concluding  with  2^  =  8  tiles  around  each  of  the  poles. 
These  spherical  wireframes  are  usable  without  modification  in  problems  dealing  with  spherical 
plumes.  Figure  5  is  the  m  =  4  case  with  a  diameter  of  1  m.  Actually,  wireframes  of  all  shapes  are 
generated  in  this  “unit”  form  (with  a  maximum  diameter  of  1  m)  for  convenience  and  are  simply 
scaled  up  or  down  as  required  for  various  problems. 

To  produce  a  wireframe  for  a  nonspherical  plume,  we  first  simply  project  the  spherical 
wireframe’s  vertices  radially.  Figure  6  is  such  a  radial  projection  of  the  wireframe  of  Fig.  5  onto  an 
ellipsoidal  surface  whose  x,  y,  and  z  dimensions  are  in  the  ratio  1:2:3.  Since  this,  too,  is  a  “unit” 
wireframe,  its  actual  measurements  are  1/3  m  by  2/3  m  by  1  m.  Projection  always  deforms  the 
wireframe  to  some  degree,  shrinking  some  tiles  and  stretching  others.  This  should  be  harmless 
unless  some  of  the  tiles  end  up  larger  than  a  fraction  of  the  interior  wavelength.  That  kind  of 
distortion  can  degrade  the  capability  for  coherent  simulation.  When  it  occurs,  we  can  compensate 
by  further  adjusting  the  wireframe.  We  do  this  by  a  process  of  simulated  annealing  applied  so 
as  to  minimize  the  normalized  variance  of  the  tile  edge  lengths  (Appendix  D).  This  encourages 
the  tiles  to  become  as  nearly  equilateral  and  uniform  in  area  as  the  shape  of  the  plume  and  the 
topology  of  the  wireframe  will  allow.  The  results  of  the  procedure  can  be  seen  in  Fig.  7. 

The  potential  for  modeling  with  such  a  wireframe  can  be  assessed  from  the  relation 


P  a(/3ci„t)2 

(obtained  in  Appendix  E)  in  which  a  and  are  the  plume’s  area  and  average  interior  sound  speed 
and  a  and  are  the  constants  0.43  and  0.25,  respectively.  For  a  given  plume,  this  expression 
pertains  to  a  uniform  equilateral  wireframe  with  tile  sides  equal  to  Aint/4.  It  can  be  used  to 
estimate  either  the  minimum  N  required  for  a  particular  /  or  the  maximum  /  at  which  a  given 
N  is  feasible.  We  will  use  it  for  both  purposes  and  will  refer  to  the  result  as  the  “marginal  tiling 
estimate”  for  either  N  or  f  according  to  context. 
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Fig.  6  -  Spherical  wireframe  of  Fig.  5  projected 
onto  an  ellipsoidal  plume  shape 


We  note  as  an  aside  that,  although  the  above  tiling  scheme  worked  flawlessly,  it  is  admittedly 
somewhat  lacking  in  flexibility  in  that  it  provides  no  intermediate  possibilities  between  =  16 
tiles  for  m  =  3  and  N  =  192  for  m  =  4.  While  the  latter  proved  entirely  satisfactory  here,  it 
does  represent  a  degree  of  numerical  “overkill”  that  is  undesirable  for  two  reasons.  First,  there 
is  the  practical  matter  of  computer  run  time.  This  increases  roughly  as  so  that  using  an 
m  =  4  wireframe,  for  example,  when  the  marginal  tiling  estimate  is  only  60  means  unnecessarily 
prolonging  the  execution  by  a  factor  of  (192/60)^  «  10.  Secondly,  a  tiling  method  that  allowed  any 
integral  N,  rather  than  just  16  or  192,  would  open  the  way  to  using  a  method  like  Richardson’s 
deferred  approach  to  the  limit  [11]  to  estimate  the  values  of  integrals  in  the  iV  — ►  oo  limit.  This 
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could  provide  increased  accuracy  for  less  effort.  In  the  end,  it  proved  prohibitively  difficult  to 
devise  such  a  scheme  that  achieved  the  desired  flexibility  while  covering  the  sphere  with  tiles  of 
controlled  size  and  shape. 


5.2  Numerical  Computation  of  the  Scattering  Amplitude 


Equation  (25)  is  a  2i\f -dimensional  linear  problem  of  the  form  Mx  =  b  and  could  be  solved  by 
a  variety  of  standard  methods.  We  chose  the  LU  decomposition  method  with  iterated  improvement 
—  a  complex  version  of  the  implementation  in  Numerical  Recipes  [11].  For  the  modeling,  we  allow  a 
discrete  set  of  sources,  i.e.,  incident  multipath  directions  {Sf^;  cr  =  1,  •  •  • ,  Nsrc},  each  corresponding 
to  a  linear  problem 

Mx,,  =  b„  ,  (31) 


where 


P(s<t)  \ 
Q{s.)  )- 


(32) 


For  a  discrete  set  of  outgoing  directions  {fp;  p  =  1,  •  •  • ,  Nrec},  Eq.  (29)  amounts  to 

F  pa  =  CpXa 

where  Fpa  =  F(fp|  -  s„)  and 

Cp  = 

When  there  are  Ngrc  >  1  coherent  arrivals,  the  overall  scattering  amplitude  has  the  form 

Narc 

Fp  —  ^  OtaFpa 
<7=1 


(33) 

(34) 


(35) 


with  the  relative  amplitudes  and  phases  of  the  multipath  signals  specified  by  the  complex  Qc 
coefficients. 


Since  the  whole  procedure  is  analytically  equivalent  to 

/^arc 

Fp  =  cjM-^  ( 

\ct=1 

numerical  difficulties  can  be  expected  if  M  becomes  ill-conditioned.  We  encountered  this  problem 
(and  were  forced  to  switch  from  LU  to  singular  value  decomposition  as  a  remedy)  only  for  environ¬ 
mental  parameter  values  far  outside  the  physical  range  for  this  application  —  interior  sound  speeds 
of  50  m/s  for  instance.  Realistic  model  plumes  did  not  resonate  with  a  high  enough  Q- factor  to 
cause  the  problem.  Instead,  the  main  problem  for  our  approach  is  that  2N  may  be  quite  large. 
This  requires  plenty  of  RAM  or  a  virtual  memory  system.  Fortunately,  current  computers  are 
adequate  for  most  practical  cases.  The  details  of  this  point  are  discussed  in  Appendix  E. 

BIRPS  is  implemented  in  standard  FORTRAN  77  and  follows  the  pseudocode  outline  below. 
Steps  (a)  through  (c)  define  the  problem  analytically  and  produce  a  discrete  representation  of  it. 
Step  (d)  performs  the  LU  decomposition.  When  Ngrc  >  L  loops  (e)  and  (1)  scatter  the  multipath 
arrivals  coherently. 
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(a)  select  the  shape  of  S  and  construct  the  wireframe 

(b)  set  the  properties  of  £,  1 

(c)  compute  M 

(d)  LU-decompose  M 

(e)  for  (T  from  1  to  Ngrc 

(f )  compute  6(7 

(g)  compute  Xcr  by  LU  back-substitution 

(h)  apply  iterated  improvement  to  =  0 

next  (T 

(i)  for  p  from  1  to  Nrec 

(j)  initialize:  Fp  =  0 

(k)  compute  Cp 

(l)  for  a  from  1  to  Ngj-c 

(m)  compute  Fpa  by  Eq.  (33) 

(n)  accumulate :  Fp  =  Fp  +  *  Fpo 

next  a 

next  p 

The  plume  shape  is  specified  in  step  (a)  by  a  ‘unit  radius  function’  subroutine 

SUBROUTINE  urf  (  theta,  phi,  urad,  uonv  ) 

whose  inputs  theta,  phi  are  the  polar  angles  of  a  point  on  the  plume  surface  and  whose  outputs 
urad,  uonv  are  the  radius  and  unit  outward  normal  vector  at  that  point.  An  interior  sound  speed 
profile  is  specified  in  terms  of  the  slowness  c~^  by  a  real- valued  function  svp  whose  entry  points 
are  the  slowness  itself, 

ENTRY  slow  (xi), 

and  its  three  spatial  gradient  components, 

ENTRY  gam  x  (xi), 

ENTRY  gam  y  (xi),  and 
ENTRY  gam  z  (xi). 

The  input  argument  xi  is  a  distance  which  specifies  the  physical  location  in  terms  of  a  reference 
point  and  direction  which  are  communicated  in  a  COMMON  block.  The  distance,  reference  point, 
and  direction  correspond  to  y,  and  u,  respectively,  in  Eq.  (20),  and  the  functions  are  evaluated 
at  the  point  z(^). 

Once  the  plume’s  shape  and  interior  sound  speed  are  specified,  the  numerical  computations 
closely  follow  the  analysis  laid  out  in  preceding  sections.  In  particular,  the  Green’s  functions  are 
computed  as  indicated  in  Section  3  with  the  interior  wavenumber  gradient  g,  Eq.  (23),  obtained 
directly  from  the  components  of  the  slowness  gradient,  and  the  integrals  e  and  f,  Eq.  (24),  com¬ 
puted  using  Romberg  quadrature  routines  on  the  components  of  g.  Thus,  to  implement  additional 
plume  shapes  and  interior  sound  speeds  beyond  those  used  in  this  report,  the  user  need  only  supply 
new  versions  of  urf  and  svp. 

6,  PLUME  RESULTS 

In  the  following  subsections  we  simulate  plume  scattering  in  a  series  of  situations.  Except  as 
noted,  all  of  these  situations  incorporate  the  pressure-relief  condition  at  the  sea  surface  and  involve 
a  single  arrival.  We  will  be  computing  the  scattering  amplitude  F  and  plotting  it  (or  rather  its 
modulus,  since  F  is  complex)  as  a  function  of  various  parameters  in  the  following  figures.  F  is 
referred  to  a  unit-amplitude  incident  plane  wave  and  has  the  physical  dimensions  of  a  length. 
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The  sound  speed  at  any  point  inside  a  plume  can  be  related  to  the  air  volume  fraction 
through  the  general  expression  [2,  Section  8.9] 


c  = 


{^Pg  +  (1  ~  0)pe) 


0  1-0 


^Pgcj 


+ 


-1/2 


(36) 


which  gives  the  effective  sound  speed  c  for  a  suspension  of  gas  bubbles  in  a  liquid  in  terms  of  /3,  the 
volume  fraction  of  gas.  Expert  opinion  [12]  currently  distinguishes  three  bubble  structure  categories 
by  air  volume  fraction,  namely  “tenuous”  {0  <  10““*),  “intermediate”  (10“'*  <  0  <  10“^),  and 
“dense”  (10“^  <  0).  The  names  refer  to  the  number  density  of  air  bubbles  thought  to  be  present. 
Figure  8  is  computed  from  Eq.  (36)  for  air  in  water  and  covers  the  full  “intermediate”  range.  The 
following  simulations  deal  with  interior  sound  speeds  in  the  300-800  m/s  range  and  hence  with 
the  class  of  bubble  structures  known  as  intermediate  plumes.  The  same  method  would  also  apply 
to  tenuous  plumes,  without  reliance  on  the  approximations  usually  invoked  in  that  regime  [13-15], 
and  to  dense  plumes,  provided  they  are  static  and  their  sound  speed  gradients  are  not  too  large. 


Fig.  8  -  Sound  speed  from  Ekj.  (36)  for  air  suspended  in  water  throughout  the  “intermediate 

plume”  range,  0.01%  <  0  <  1% 


6.1  Fundamental  Features  in  Simple  Situations 

We  begin  by  applying  the  BIRPS  model  to  simple  situations  involving  a  spherical  plume 
with  an  isovelocity  interior.  An  important  objective  of  this  first  set  of  simulations  is  to  validate 
the  model  using  fairly  realistic  physical  parameters. 

6.1.1  Azimuth  Dependence 

Figure  9  examines  two  things:  the  azimuth  dependence  of  the  scattering  and  the  number 
of  tiles  required  for  accurate  computations.  The  incident  azimuth  is  22.5°  (i.e.,  360°  16)  so 
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Fig.  9  -  Scattering  amplitude  eis  a  function  of  azimuth  at  125  Hz 

that  both  the  m  =  3  wireframe  (with  eight  equatorial  tiles)  and  the  m  =  4  wireframe  (with  16) 
will  be  symmetric  in  azimuth  about  the  source  direction.  Thus  the  scattering  amplitude  as  a 
function  of  receiver  azimuth  should  inherit  that  same  symmetry.  The  figure  does  indeed  have  the 
expected  symmetry  relative  to  the  source  direction  (the  broken  straight  line).  Given  the  diameter, 
frequency,  and  interior  sound  speed,  the  marginal  tiling  estimate  (from  Eq.  (30),  see  Appendix 
E)  for  the  minimum  number  of  tiles  required  is  29.  Indeed,  the  m  =  3  case  (dashed),  which  has 
only  16  tiles,  is  about  10%  below  the  m  =  4  case  (solid)  with  192.  The  scattering  amplitude  is 
small  (only  about  0.04)  for  these  near-grazing  source  and  receiver  elevations  (5°  and  2°  below  the 
horizontal,  respectively).  We  will  see  that  this  is  because  these  elevations  lie  in  a  horizontal  notch 
in  the  vertical  dipole-type  beam  pattern  formed  by  the  plume  and  its  surface  image. 
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6.1.2  Resonance  Spectrum 

Figure  10  examines  the  frequency  dependence  of  scattering  from  the  same  plume  without 
the  pressure-release  interface.  Here  the  relevant  angle  is  simply  the  total  scattering  angle  x)  = 
arceos(— s-r),  which  is  taken  to  be  45°.  Results  from  the  present  computer  model  are  shown  as 
circles  —  o  for  m  =  3  and  •  for  m  =  4.  The  heavy  line  represents  the  scattering  amplitude  as 
computed  by  the  semi-analytic  partial  wave  method  described  in  Appendix  F.  The  broken  lines  are 
the  individual  ^  =  0,  ■  •  • ,  3  partial  wave  contributions  —  s,  p,  d,  f-wave  scattering  in  atomic  physics 
jargon.  The  peaks  clearly  correspond  to  scattering  resonances.  The  ^th  one  lies  very  close  to  the 
frequency  (£  -I-  1)  x  125  Hz  that  would  be  anticipated  from  a  “back  of  the  envelope”  calculation 
using  Celt  =  oo  and  it  can  be  expected  to  exhibit  the  familiar  {£  -I-  l)-lobed  structure  in  the  angle 
6  (Appendix  F).  The  Q-factor  (i.e.,  the  peak  frequency  divided  by  the  half- width)  increases  with 
frequency:  Q  «  1, 3, 10  for  ^  =  0, 1, 2. 


FREQUE(«:y  (Hz) 


Fig.  10  -  Model  results  for  m  =  3  (o)  and  m  =  4  (•).  For  compeurison,  computations  from 
Appendix  F;  |F|  (heavy  line)  and  [Ftl  for  £  =  0, 1,2,3  (broken  lines)  as  functions  of  / 


A  marginal  tiling  estimate  indicates  that  the  m  =  3, 4  contributions  should  be  accurate  up 
to  92  Hz  and  320  Hz,  respectively.  According  to  the  figure,  however,  this  kind  of  estimate  is 
overly  conservative.  The  o  symbols  depart  from  the  heavy  line  somewhere  beyond  195  Hz  and  the 
•  symbols  are  adequately  close  to  it  through  at  least  405  Hz.  The  numerical  model  exceeds  the 
marginal  tiling  expectations,  failing  only  with  the  the  onset  of  the  ^  =  3  resonance  around  445  Hz. 
This  can  be  understood  by  considering  the  radiation  pattern  for  this  resonance.  It  has  four  lobes 
in  the  range  0  <  <  tt  while  the  wireframe  (Fig.  5)  has  only  eight  “latitude  bands”  of  tiles  in  that 

interval.  That  angular  undersampling  —  only  two  tiles  per  angular  lobe  —  appears  to  cause  the 
breakdown  at  this  point.  The  method  should  work  at  least  this  well  for  smaller  or  faster  plumes, 
since  their  resonances  will  be  shifted  to  higher  frequencies  [16,  Chapter  VI,  Section  2]. 


18 


NRL  REPORT  9391 


6.1.3  Elevation  Dependence 

In  Fig.  11  we  examine  the  scattering  response,  at  0°  elevation  (downward)  and  45°  elevation 
in  the  forward  azimuth  plane,  for  the  Fig.  9  setup  in  the  frequency  band  105  to  145  Hz.  The 
s-wave  resonance  is  clearly  present  here,  too,  although  comparison  of  Figs.  9  and  11  to  Fig.  10 
indicates  that  more  tiles  are  needed  to  model  it  when  the  pressure-release  interface  is  present.  In 
recognition  of  this,  the  remainder  of  the  modeling  will  emphasize  m  =  4  wireframes. 
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Fig.  11  -  Scattering  amplitude  a.s  a  function  of  frequency  near  125  Hz  for  both  0°  and  45° 

receiver  elevation 
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For  the  same  situation,  Fig.  12  shows  the  dependence  on  output  elevation  in  the  entire 
vertical  half-plane  containing  the  source  direction.  The  source  direction  is  indicated  as  usual  by 
the  broken  line.  Back-scattering  corresponds  to  the  right  half  of  the  figure  and  forward  scattering 
is  on  the  left.  The  dashed  curves  are  for  frequencies  above  the  resonance  peak.  All  of  the  curves 
have  a  high  degree  of  forward/backward  symmetry,  evidently  because  onlj'  the  symmetric  s-wave 
resonance  contributes  in  this  frequency  range. 
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Fig.  12  -  Scattering  amplitude  as  a  function  of  elevation,  {parameterized  by  frequency.  The 
discrete  frequencies  are  the  same  as  in  Fig.  11:  105  Hz  to  145  Hz  in  5  Hz  steps.  Solid: 
/  <  130  Hz,  dashed:  130  Hz  <  /. 
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6.1.4  Depth  Dependence. 

For  the  same  plume,  Fig.  14  shows  the  near-vertical  scattering  response  to  a  near-vertical 
arrival  as  a  function  of  the  plume’s  center  depth  d  in  the  interval  2  m<  d  <4  m.  For  comparison, 
consider  what  Fig.  14  would  look  like  if,  rather  than  this  large  resonant  plume,  we  were  dealing  with 
a  small  noiiresonant  object  —  a  “point  target.”  Then,  because  both  elevations  are  nearly  vertical, 
the  d  dependence  would  be  sixi^ (kextd)  ■  This  vanishes  at  the  surface  and  at  d  =  =  6  m  and 

is  symmetric  about  its  peak  at  d  =  3  m.  The  actual  curve  in  Fig.  14  is  similar  to  that  in  functional 
form  but  has  its  peak  shifted  upward  by  about  0.75  m,  presumably  due  to  the  plume’s  finite  size 
and  internal  structure. 
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Next  we  turn  to  the  near-grazing  geometry  where,  in  view  of  its  greater  practical  importance, 
we  go  into  more  detail.  The  source  and  receiver  are  held  at  grazing  angles  of  2.5°  and,  except 
as  noted,  other  parameters  remain  as  indicated  in  the  caption  to  Fig.  14.  The  object  is  again 
to  determine  the  depth  dependence  of  the  response  in  the  interval  2  m  <  d  <  4  m.  This  was 
done  through  BIRPS  simulations  at  0.25  m  intervals.  First,  as  a  benchmark  case,  the  plume 
diameter  D  was  reduced  to  only  0.5  m.  The  plume  is  small  relative  to  the  exterior  wavelength 
{D  -C  Xeit)  SO  that,  as  d  — *  4  m,  the  response  should  exhibit  the  quadratic  depth  dependence 
characteristic  of  a  point  target  in  this  grazing  geometry.  The  simulations  accurately  confirm  the 
expected  S  dependence.  For  the  second  simulation  series,  the  2  m  diameter  was  reinstated  and 
the  result  was  found  to  be  more  nearly  linear  in  d.  Since  this  plume  w^as  resonant  at  the  frequency 
considered  here,  it  was  not  clear  whether  this  nonquadratic  depth  dependence  was  produced  by 
the  resonance  or  simply  by  the  larger  size.  The  third  simulation  series  addressed  that  question. 
Here,  the  interior  sound  speed  was  increased  to  1000  m/s,  which  shifted  the  lowest  resonance  well 
above  the  operating  frequency.  The  result  for  this  large  nonresonant  plume  was  roughly  linear 
for  d  ss  2  m  but  became  quadratic  for  d  >  4  m.  Figure  15  summarizes  all  three  near-grazing 
simulations.  If  it  is  assumed  that  the  depth  dependence  follows  some  power  law  |F(d)|  oc  (F,  then 
the  exponent  can  be  assessed  from  the  simulations  by  plotting  p  =  |F(d)|'d/|F(d)|  vs  d.  This  is 
done  in  Fig.  15  with  the  depth  derivative  estimated  by  finite  differences.  The  small-plume  result 
is  the  top  curve  which,  as  expected,  indicates  p  =  2.  The  middle  curve  corresponds  to  the  large 
nonresonant  plume.  Its  depth  dependence  starts  out  as  rather  linear  at  d  =  2.25  m  but  becomes 
effectively  quadratic  by  d  =  3.75  m.  The  bottom  curve  is  the  result  for  the  large  resonant  plume. 
The  behavior  is  more  nearly  linear  for  d  <  4  m,  although  the  trend,  if  extrapolated,  is  toward 
quadratic  behavior  at  greater  depths.  Thus,  the  deviation  from  quadratic  dependence  seems  to 
be  in  part  a  result  of  large  D/d.  Resonance,  however,  enhances  it,  possibly  by  giving  the  plume 
a  larger  effective  diameter.  Of  course,  a  finite  diameter  and  a  resonant  interior  do  more  than  just 
alter  the  depth  dependence  of  a  plume’s  scattering  response.  They  will  both  raise  the  absolute 
level  too.  as  shown  in  Fig.  16.  The  large  nonresonant  plume  scatters  at  levels  10-12  dB  greater 
than  the  '‘point-target”  plume  and  resonance  adds  another  10-12  dB. 


6.1.5  Summary 

We  find  that  the  low-frequency  scattering  responses  of  spherical  model  intermediate  plumes 
contain  multiple  weak-to-moderate  resonances.  The  lowest-frequency  resonance  is  an  s-wave  fea¬ 
ture  with  a  uniform  azimuth  response  and  an  elevation  response  similar  to  that  of  a  vertical  dipole. 
The  next  resonance  has  a  p-wave  character  so  that  its  azimuth  response  has  forward  and  backward 
lobes  and  its  elevation  response  is  quadrupole-like.  As  the  frequency  increases  through  this  reso¬ 
nance,  the  nature  of  the  response  shifts  from  weakly  back  scattering  to  strongly  forward  scattering. 
Thus,  for  a  receiver  in  the  forward  direction  and  an  incident  signal  with  a  power  spectrum  concen¬ 
trated  around  the  second  resonance  frequency,  the  scattered  signal  would  be  high-pass  filtered.  It 
is  also  clear  that  the  depth  dependence  of  .scattering  from  realistic  intermediate  plumes  in  the  first 
few  meters  below  the  surface  is  very  different  from  what  would  be  expected  from  a  point  target 
-  especially  at  resonance  frequencies. 

6,2  Effects  of  Multiple  Arrivals 

Now  we  turn  to  an  application  in  .shallow  water  with  .several  propagation  multipaths.  For 
this  simulation,  we  used  eight  vertical  arrivals  from  a  source  in  a  shallow  Pekeris  environment. 
Figure  17  summarizes  the  simulation  parameters  and  the  eight  multipaths.  These  multipaths  are 
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Fig.  15  -  Power-law  analysis  of  the  depth  dependence  of  the  scattering  amplitude  for  near¬ 
grazing  elevations  at  125  Hz 
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Fig.  16  -  Scattering  amplitude  for  near-grazing  elevations  as  a  function  of  plume  center  depth 
at  125  Hz.  Vertical  axis  is  in  dB  relative  to  the  value  for  D  =  0.5  m,  d  =  2  m. 
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Fig.  17  -  Eigenray  multii>aths  and  simulation  parameters 


the  eigenrays  that  undergo  total  reflection  at  the  sediment  in  this  geometry.  The  partially  reflected 
rays  all  had  power  levels  at  least  16  dB  lower  and  have  been  neglected.  Figure  18  shows  the  azimuth 
dependence  of  the  scattering  at  the  lowest  resonance  frequency.  This  figure  is  a  multipath  version 
of  Fig.  9.  It  has  the  necessary  symmetry  but  the  amplitude  is  larger  in  this  case  by  a  factor  of 
more  than  twenty,  which  may  at  first  glance  seem  unreasonably  high.  Indeed,  even  if  all  of  the 
rays  somehow  managed  to  arrive  at  the  direct  path  elevation  (87.71°)  and  combine  in  phase,  that 
could  accoimt  for  no  more  than  a  factor  of  eight.  It  turns  out  that  most  of  this  increase  is  caused 
by  the  strong  dependence  of  the  scattering  amplitude  on  source  elevation.  The  rays  that  produce 
the  greatest  response  are  those  with  the  steepest  grazing  angles.  In  f^lct,  separate  simulations  have 
shown  that  the  steepest  ray  (elevation  59.04°  and  4  bottom  bounces)  accounts  for  60%  of  the  total 
scattering  response  in  Fig.  18  while  the  near-grazing  direct  ai-rival  (87.71°  elevation)  produces  only 
3%.  It  may  be  worth  noting  here  that,  for  shorter  ranges,  even  some  of  the  partially  reflected  rays 
may  need  to  be  included. 

Figure  19  shows  the  scattering  response  as  a  function  of  output  azimuth  at  the  second 
resonance  frequency.  Here,  as  in  Fig.  13,  there  are  distinct  forward  and  backward  lobes  due  to 
the  ‘p-wave’  nature  of  the  resonance  and  the  proximity  of  the  pressure-release  boundary.  This  is 
further  illustrated  in  Figs.  20  and  21  which  show  the  elevation  dependence  at  the  first  and  second 
resonances. 

6.3  Effects  of  Nonspherical  Plume  Shape 

In  this  section  we  consider  the  scattering  response  of  the  ellipsoidal  plume  shown  in  Fig.  6. 
The  semiaxes  are  aligned  with  the  standard  x,  y,  z  directions,  but  their  lengths  are  in  the  pro¬ 
portions  1:2:3  so  that  this  is  not  a  figure  of  revolution.  As  before,  the  largest  overall  dimension 
is  2.0  m,  so  we  are  dealing  with  a  plume  that  stands  2  m  tall  (in  the  z  direction),  has  a  width  (y 
direction)  of  Ij  m  and  is  |  m  thick  (x  direction). 
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Fig.  18  -  Scattering  amplitude  as  a  function  of  azimuth  at  125  Hz.  This  is  a  version  of 
Fig.  9  with  eight  arrivals.  For  simplicity,  only  the  m  =  4  result  is  plotted. 
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Fig.  19  -  Scattering  amplitude  as  a  function  of  azimuth  at  240  Hz 
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Fig.  20  -  Scattering  amplitude  as  a  function  of  elevation  at  125  Hz 
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Figure  22  examines  the  resonance  spectrum  for  scattering  from  this  structure.  To  avoid 
unnecessary  complications,  this  is  done  without  the  pressure-relief  interface.  The  model  results  are 
plotted  as  circles  —  o  for  the  nonannealed  wireframe  of  Fig.  6  and  •  for  the  annealed  wireframe 
of  Fig.  7.  Since  this  plume  is  physically  smaller  and  less  symmetric  than  the  preceding  spherical 
one,  its  resonant  eigen-frequencies  should  be  higher  and  less  degenerate  [16,  Chapter  VI,  Section 
2].  This  expectation  is  supported  by  the  vertical  dashed  lines  that  estimate  the  lowest  few  reso¬ 
nance  frequencies  according  to  analytic  formulas  from  Morse  and  Feshbach  [17,  page  1420].  These 


formulas. 

/~ -^(0.7655 +  0.191e)  • 

•  •  B/A  «  1 

(37) 

and 

••  B/A  <  1, 

(38) 

yield  the  frequencies  listed  in  Table  1.  They  actually  represent  estimates  for  the  lowest  resonance 
of  a  2-D  elliptical  region  with  semi-axes  A,B{B<  A)  and  eccentricity  e  =  y^l  —  (B/A)'^.  For  the 
y-z  plane,  for  example,  BIA  =  2/^  =  0.67.  The  interior  sound  speed  is  c  and  the  exterior  sound 
speed  is  effectively  infinite.  Equation  (37)  applies  to  a  nearly  circular  ellipse  and  Eq.  (38)  to  a 
highly  eccentric  one.  Two  observations  may  be  made  from  Fig.  22.  First,  annealing  is  important 
since,  without  it,  the  result  computed  for  the  lowest  resonance  is  about  20  %  too  high.  Secondly, 
the  locations  of  the  resonances  agree  with  the  predictions  from  Morse  and  Feshbach  [17]  about  as 
well  as  could  be  expected,  considering  that  the  latter  derive  from  a  2-D  depiction  of  the  problem 
with  Celt  =  oo-  The  lowest  resonance,  for  example,  is  about  20  Hz  below  the  256  Hz  estimate  from 
Eq.  (38).  The  nature  of  the  peak  at  540  Hz  is  not  absolutely  certain.  A  marginal  tiling  estimate 
indicates  that  the  modeling  results  for  this  wireframe  should  be  reliable  up  to  at  least  475  Hz. 
However,  since  this  type  of  estimate  has  proved  too  conservative  in  the  spherical  case,  the  peak 
in  question  may  well  be  an  authentic  spectral  feature  analogous  to  a  p-wave  peak  for  spherical 
plumes;  i.e.,  a  genuine  resonance  which  is  simply  not  predicted  by  Eqs.  (37-38)  because  it  is  of 
higher  order. 

Table  1  —  Approximate  Resonance  Frequencies 
from  Morse  and  Feshbach 


Plane  of 
Ellipse 

B/A 

Frequer 
From  Eq.  (37) 

cy  (Hz) 

From  Eq.  (38) 

y-z 

0.67 

454 

256 

x-z 

0.33 

473 

486 

x-y 

0.50 

698 

497 

Figures  23  and  24  use  the  annealed  wireframe  and  incorporate  the  pressure-relief  interface. 
Like  Figs.  9  and  12,  they  are  done  at  the  lowest  resonance  frequency  of  the  plume  concerned. 
However,  because  this  plume  is  not  symmetric  with  respect  to  the  incident  azimuth.  Fig.  23  need 
not  retain  the  circular  symmetry  of  Fig.  9.  A  close  inspection  bears  this  out.  The  figure  is  fiattened 
by  a  factor  of  0.94  in  y  direction  and  exhibits  reflection  symmetry  relative  to  the  x  and  y  axes 
rather  than  to  the  source  direction.  This  is  normal  for  the  response  of  an  object  driven  at  an 
isolated  resonance.  Although  the  power  level  of  the  resonant  radiation  will  certainly  depend  on 
the  direction  from  which  the  object  is  excited,  the  directional  pattern  will  depend  only  on  the 
structure  of  the  object.  Figure  23  inherits  its  x-y  symmetry  from  that  of  the  plume.  Likewise,  the 
elevation  display  in  Fig.  24  retains  the  forward/backward  symmetry  of  the  plume  itself.  Figure  24 
also  confirms  the  impression  made  by  Fig.  22  that  this  lowest  resonance  is  a  lower  Q  feature  for 
the  flattened  ellipsoid  than  for  the  sphere. 
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Fig.  22  -  Scattering  amplitude  as  a  function  of  frequency.  Dashed  lines  indicate  analytic 
estimates  of  the  resonance  frequencies  as  described  in  the  text. 
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Fig.  24  -  Scattering  amplitude  as  a  function  of  elevation,  parameterized  by  frequency.  The 
discrete  frequencies  are  180,  210,  240  Hz  (solid)  and  250,  270  Hz  (dashed). 
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6.4  Effects  of  the  Gross  Features  of  the  Interior  SVP 

We  continue  with  the  ellipsoidal  plume  envelope  from  the  previous  section,  but  here  we  allow 
two  features  of  the  interior  sound  speed  —  the  central  value  cq  and  the  vertical  gradient  o'  —  to 
be  adjustable  parameters. 

6.4-1  Weak  Vertical  Gradient 

The  gradient  will  be  a  constant  c'  =  50  (m/s)/m.  Thus,  relative  to  the  sound  speed  cq 
at  the  center  of  the  plume,  the  medium  is  50  m/s  slower  at  the  top  and  50  m/s  faster  at  the 
bottom.  Without  firm  data  about  the  gradients  actually  present  in  intermediate  plumes  in  the 
oceein,  we  are  unable  to  be  any  more  realistic  about  it.  However,  since  buoyancy  should  concentrate 
more  air  at  the  top,  this  is  at  least  a  step  in  the  right  direction.  For  this  case,  Eq.  (18)  yields 
£wkb  ^  d /oj  ~  0.03  so  that  the  use  of  Eq.  (19a)  for  the  interior  Green’s  function  is  well  justified. 

The  first  computation  uses  a  sound  speed  of  cq  =  500  m/s  at  the  center.  Figure  25  examines 
the  impact  on  the  azimuth  dependence  of  the  scattering  at  240  Hz.  The  solid  curve  is  produced 
with  the  annealed  wireframe  of  Fig.  7.  It  is  indistinguishable  from  the  result  in  Fig.  23  which  has 
the  same  500  m/s  sound  speed  everywhere  without  any  gradient.  If  a  gradient  of  this  size  has  any 
effect  at  all,  it  is  not  to  be  found  at  near-grazing  elevations.  It  was  anticipated  that  annealing 
might  b;  ^  especially  important  here  because  its  effect  is  to  reconfigure  the  tiles  of  the  upper  half  of 
the  wireframe  so  as  to  improve  coverage  near  the  top  where  the  shortest  interior  wavelengths  occur. 
The  dashed  curve  was  computed  using  the  unannealed  wireframe  of  Fig.  6.  Evidently,  annealing 
is  important  here  too,  although  no  more  so  than  in  the  const£mt-svp  case  of  Fig.  22.  Figure  26 
shows  the  effect  of  the  sound  speed  gradient  on  the  elevation  and  frequency  dependence.  Again, 
the  sound  speed  gradient  makes  virtually  no  difference.  There  are  only  a  few  minor  deviations 
from  Fig.  24,  mainly  at  higher  frequencies.  The  conclusion  is  clearly  that  interior  sound  speed 
gradients  in  this  regime  are  unimportant. 

6.4.2  Reduced  Central  Speed 

Figure  27  examines  the  impact  of  the  centred  sound  speed  on  the  azimuth  dependence  — 
again  at  240  Hz.  The  weak  sound  speed  gradient  is  retained  but  the  speed  at  the  center  of  the 
plume  is  reduced  to  cq  =  400  m/s.  Again,  the  solid  curve  is  produced  with  the  annealed  wireframe. 
This  curve  differs  from  the  one  in  Fig.  25  in  two  important  respects:  (a)  the  level  of  the  scattered 
radiation  is  reduced  and  (b)  although  the  entire  plume  is  still  reflection-symmetric  with  respect 
to  the  X  and  y  axes,  its  scattering  response  no  longer  is.  Both  phenomena  can  be  explained 
by  the  fact  that,  due  to  its  altered  internal  structure,  this  plume  is  no  longer  being  driven  at 
resonance.  Whenever  a  driven  system  is  detuned  firom  a  resonance,  its  output  level  is  lowered  by 
an  amount  related  to  the  Q  of  the  resonance  and  the  amount  of  detuning.  The  asymmetry  here  is 
to  be  expected  for  any  case  that  is  not  dominated  by  a  single  resonance.  When  the  frequency  lies 
between  two  resonance  peaks,  they  both  affect  the  scattering  response  and  neither  one  can  fully 
impose  its  characteristic  symmetry  on  the  result. 

Figure  28  shows  the  effect  of  the  lower  central  sound  speed  on  the  elevation  and  frequency 
dependence.  The  “groundstate”  resonance,  which  was  in  the  210  to  240  Hz  range  when  cq  was 
500  m/s,  has  been  shifted  down  to  at  least  180  Hz  —  roughly  in  proportion  to  co«  as  would  be 
expected.  As  anticipated  from  its  slight  asymmetry,  the  240  Hz  azimuth  response  shown  in  Fig.  27 
is  somewhat  off  resonance  (by  roughly  30  to  60  Hz). 
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Fig.  25  -  Scattering  amplitude  as  a  function  of  azimuth  at  240  Hz.  This  is  a  version  of 
Fig.  23  with  a  vertical  gradient  in  the  interior  svp,  both  with  and  without  annealing  in 
the  wireframe  construction. 
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Fig.  26  -  Scattering  amplitude  as  a  function  of  elevation,  parameterized  by  frequency.  The 
discrete  frequencies  are  180,  210  and  240  Hz  (solid)  and  250  and  270  Hz  (dashed).  This  is 
similar  to  Fig.  24  but  with  a  we^lk  vertical  gradient  in  the  interior  svp. 
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Fig.  27  -  Scattering  amplitude  a.s  a  function  of  azimuth  at  240  Hz.  This  is  a  version  of 
Fig.  2.5  with  a  lower  central  sound  speed,  both  with  and  without  annealing  in  the  wireframe 
construction. 


Fig.  28  -  Scattering  amplitude  as  a  function  of  elevation,  parameterized  by  frequency. 
The  discrete  frequencies  are  150  and  180  Hz  (solid)  and  190  and  240  Hz  (dashed).  This  is 
similar  to  Fig.  26  with  a  lower  central  sound  speed. 


6.4-3  Strong  Vertical  Gradient 


We  now  return  cq  to  its  original  value  of  500  m/s  and  apply  the  strongest  gradient  that  we  can 
accommodate  with  confidence  in  the  neighborhood  of  240  Hz.  First  we  adopt  a  maximum  allowable 
value  of  0.2  for  the  parameter  ewKB  ^  to  be  used  in  the  modeling.  Since  ewkb  we  are 
justified  in  using  Eq.  (19a)  for  the  interior  Green’s  function.  Then  we  use  the  msocimum  gradient 
consistent  with  that  choice  at  240  Hz;  i.e.,  d  =  0.2  x  27r  x  240  «  300  (m/s)/m.  It  may  be  seen 
in  Fig.  29  that  the  effect  of  this  strong  vertical  gradient  is  to  sharpen  the  resonance  (raise  the  Q) 
but  not  to  shift  it. 
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Fig.  29  -  Scattering  amplitude  as  a  function  of  elevation,  parameterized  by  frequency.  The 
discrete  frequencies  are  180  210  and  240  Hz  (solid)  and  250  and  270  Hz  (dashed).  This  is 
similar  to  Fig.  24  but  with  a  strong  vertical  gradient  in  the  interior  svp. 


7.  SUMMARY  AND  CONCLUSIONS 

We  have  developed  a  numerical  implementation  of  the  boundary  integral  treatment  of  scat¬ 
tering  from  acoustically  penetrable  bubble  structures  of  very  general  type.  The  underlying  theory 
includes  the  effects  of  collective  bubble  resonances  and  our  implementation  provides  options  for  the 
presence  of  a  pressure-release  sea  surface  and  for  coherent  multipath  arrivals.  The  implementation, 
in  the  form  of  the  BIRPS  software,  provides  a  numerical  modeling  component  for  the  investigation 
of  low-frequency  surface  reverberation  and  has  the  potential  for  wider  application  to  the  study  of 
resonant  scattering  from  other  complex  bodies. 
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We  have  been  able  to  compute  the  complex  scattering  amplitudes  of  model  structures  at 
frequencies  up  to  a  few  hundred  Hz  for  a  range  of  size,  shape,  depth,  and  interior  sound  speed 
profile.  These  computations  were  done  with  parameter  values  chosen  to  simulate  plumes  of  the 
“intermediate”  type.  We  find  that  such  plumes  have  multiple  weak-to-moderate  resonances  in 
the  sub-kHz  band  of  interest.  The  nonresonant  scattering  response  for  near-grazing  geometries  is 
about  10  dB  higher  than  point-target  levels  and  tends  be  linear  in  the  depth  near  2  m,  changing  to 
quadratic  by  4  m.  Resonance  adds  approximately  another  10  dB  to  the  level  and  extends  the  linear 
dependence  to  greater  depths.  Typically,  the  lowest-frequency  resonance  is  a  low-Q  feature  whose 
response  has  a  weak  frequency  dependence  with  a  vertical  dipole-like  dependence  on  elevation 
angle.  As  frequency  is  increased,  the  next  resonance  encountered  has  a  higher  Q  and  exhibits 
forward  and  backward  scattering  lobes.  As  the  driving  frequency  is  swept  upward  through  such  a 
resonance,  the  dominant  component  of  the  response  shifts  from  backscattering  to  forward  scattering 
at  a  rate  dependent  on  the  Q.  Features  with  a  significant  impact  on  the  resonance  spectrum  of  a 
plume,  and  on  the  scattering  response  in  general,  are  the  “primary  parameters”  (size,  shape,  and 
mean  interior  sound  speed)  and  the  depth.  Strong  sound  speed  gradients  enhance  the  Q  of  the 
resonant  responses,  while  weak  gradients  have  no  appreciable  effect. 

This  work  has  focused  on  the  resonance  effects  of  bubble  structures.  We  have  neglected 
“tenuous”  bubble  clouds  (the  large,  persistent  features  sustained  by  Langmuir  circulation)  for 
that  reason.  These  may  be  of  considerable  potential  importance  in  ocean  acoustics,  but  they 
are  weak,  nonresonant  scatterers.  This  neglect,  however,  was  purely  a  matter  of  emphasis.  The 
BIRPS  software  is  certainly  also  applicable  to  tenuous  clouds  and  may  be  useful  in  that  regime, 
particularly  as  a  means  of  validating  the  approximations  made  in  other  approaches. 

Two  limitations  of  the  present  technique  are  that  it  treats  the  sea  surface  as  flat  and  the 
plume  as  static.  Although  surface  curvature  on  a  scale  comparable  to  plume  dimensions  would 
need  separate  treatment,  curvature  on  larger  and  smaller  scales  could  be  addressed  with  only 
moderate  additional  effort.  For  example,  swell  might  be  simulated  as  a  time-dependent  surface  tilt 
and  small-scale  roughness  could  be  incorporated  in  a  surface  reflection  coefficient.  As  it  stands,  our 
approach  could  deal  with  scattering  from  plumes  of  the  “dense”  class  —  provided  they  remained 
static.  (Very  strong  interior  sound  speed  gradients,  if  present  in  such  plumes,  would  require  an 
improved  formulation  of  the  interior  Green’s  function,  but  the  approach  is  still  sound.)  Actual 
dense  plumes,  however,  have  such  short  lifetimes  that  the  present  CW  implementation  of  BIRPS 
is  probably  not  adequate  for  them.  Such  fast-evolving  structures  would  require  a  broadband 
implementation  to  represent  their  time-dependent  scattering.  Incorporating  surface  swell  and 
roughness  and  broadband  signal  capability  into  BIRPS  are  promising  areas  for  future  work. 

The  work  reported  in  this  publication  was  not  done  in  a  vacuum.  Numerous  sources  have 
proved  valuable  to  the  authors  in  acquiring  the  background  to  tackle  the  subject  and  a  few  have 
provided  inspiration  for  the  particular  approach  adopted  here.  We  have  included  a  complete 
account  of  this  in  Appendix  G  to  help  put  our  work  in  context  and  to  assure  that  appropriate 
credit  is  given. 
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Appendix  A 


UNIFORM  1-D  OBJECT 


This  appendix  reviews  the  solution  to  the  one-dimensional  problem  of  CW  scattering  by  a 
uniform  object,  as  shown  in  Fig.  Al.  For  simplicity,  the  density  is  the  same  inside  and  outside  the 
scatterer.  For  a  unit-amplitude  field  incident  from  the  left,  the  solution  is  found  via  the  following 
steps: 


1.  Form  the  total  field  in  each  region  — 


field  region 

jg+ifeii  ^  <  X  <  0 

^g+tfeoa:  ^  ^g-tfcox  0  <  X  < 

j-g+tfcji  w  <  X  <  4-00  .  (Al) 

R  and  T  are  the  complex  reflection  and  transmission  coefficients. 

2.  Impose  continuity  on  the  field  and  its  x-derivative  at  x  =  0,  w. 

3.  Solve  the  resulting  equations  for  T,  R. 


km(0^C 

A 


0 


H - X 

w 


Fig.  Al  -  Sketch  of  a  1-D  problem,  showing  the  notation  used 
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It  is  easily  seen  that,  with  /x  =  6+*''°“',  the  conditions  that  accomplish  step  2  are 


we  have 


TD  = 


fcoM  -koH* 
1  1 
ko  -ko 


With  the  determinant 


M* 

feo/i 

-kon* 

1 

1 

ko 

—ko 

0 

M* 

0  - 

-fco/x* 

1 

1 

ko 

—ko 

M 

0  0 

kon 

0  0 

1 

1  - 

ko 

fci  k 

kofi 

-koM* 

1 

1 

ko 

-ko 

kon 

-kon* 

1 

1 

ko 

-ko 

0  -1 


0  -1 


-1  0 


0  -1 


1  0 


0  0 


-1  1 


=  (A:i  +  ko)^fj,* 


kofu, 


=  1li*k\{k\  -h  Atq) 


=  -2/xfci(fci  -  fco) 


=  -2i{kx  -  ko){ki  +  A:o)sin(A:ot;) 


=  Akoki 


Resonances  may  be  located  graphically  by  plotting  |i?|  or  \T\  vs  frequency.  They  should  occur 
at  approximately  the  frequencies  where  the  bound  states  occur  in  the  c  — >  oo  limit,  namely  at 
/  =  n  A/  for  n  >  0,  with  an  even  spacing  A/  =  co/2«;.  This  is  borne  out  in  the  body  of  the  report 
in  Fig.  1.  Prom  that  figure  it  is  also  clear  that,  even  without  attenuation,  the  acoustic  resonances 
have  finite  widths. 
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SURFACE  TILE  INTEGRALS 

This  appendix  evaluates  and  Sn  —  the  integrals  that  appe^l^  in  the  matrix  elements  of 
the  wireframe  continuity  condition  of  Section  4. 

Let  Vi,  V2,  V3  denote  vectors  from  the  center  of  the  plume  to  the  vertices  of  the  nth  tile. 
The  endpoint  of  the  vector  v  =  (vj  4-  V2  +  V3)/3,  then,  lies  at  the  geometric  center  of  the  tile.  We 
will  scale  v  up  or  down  as  required  so  that  its  endpoint  lies  on  the  plume  itself  and  denote  the 
resulting  vector  by  a.  Thus,  the  endpoints  of  vi,  V2,  V3  form  the  base  of  a  pyramid  whose  apex 
is  at  a  as  in  Fig.  Bl.  The  vectors  pj  =  vj  —  a.  locate  the  corners  of  that  p3Tamid  relative  to  the 
apex. 


Fig.  Bl  -  Pyramid  having  a  tile  as  its  base  eind 
its  apex  on  the  bubble  plume 

We  first  consider  the  shaded  triangular  part  of  one  face  of  the  pyramid  shown  in  Fig.  B2, 
with  r  =  (r,  0)  denoting  a  radius  vector  from  the  apex  to  a  point  lying  on  that  face,  and  compute 


M 


=  //rrfrde(l) 


over  this  region  by  adding  the  contributions  from  the  subregions  I  and  II:  M  =  Mi  +  M//  where 

re 


Ml  =  f  d6  f  dr  =  Qq 
Jo  Jo 

fp  r© 

M//  =  I  dr 

Jq  Jcos-*{q/r) 


d0 


—  {P  ~  q}^  ~  f  dr  cos  ^(q/r)  . 
Jq 


(Bl) 


(B2) 
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APEX 


Addition  €tnd  a  little  algebra  result  in 


M  =  pH{q/p) 

(B3) 

where 

(B4) 

Wn  is  the  sum  of  the  contributions  from  all  six  such  regions  (two  for  each  face). 

In  terms  of  the 

perpendicular  distance, 

„  _  IPi  X  Pj\ 

IPi-Pjl 

(B5) 

from  the  apex  to  base  edge  i-j  as  shown  in  Fig.  B3, 

(B6) 

The  roughest  estimate  for  5„  would  be  simply  the  area  of  the  tile;  i.e., 

the  area  of  the 

pyramid’s  base 

Sn  =  |piXP2  +  P2Xp3  +  P3XPll/2  . 

(B7) 

We  use  a  more  accurate  estimate. 

Sn  =  (|PlXp2|  +  |P2Xp3|  +  )P3Xpi))/2  , 

(B8) 

the  sum  of  the  areas  of  the  three  faces.  The  improvement  that  this  provides  over  Eq.  (B7)  is 
significant  when  the  height  of  the  pyramid  is  a  significant  fraction  of  its  base;  i.e.,  when  the 
number  of  tiles  is  small. 


46 


NRL  report  9391 


Fig,  B3  -  Pyramid  showing  Pj  and  qa 


for  all  three  faces 
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DIAGONAL  MATRIX  ELEMENTS 


To  obtain  the  diagonal  matrix  elements  given  in  Eq.  (28),  it  is  necessary  to  evaluate  Eqs.  (17) 
and  (19)  in  the  limit  x  Xm  where  Xm  is  the  coordinate  at  the  center  of  the  mth  tile.  In  this  limit, 
singularities  occur  in  the  integrands.  These  are  dealt  with  as  follows.  The  integrand  is  expanded 
in  a  Laurent  series  about  the  point  \x  —  Xm\  =0.  Cancellations  occur,  and  at  worst  only  a  simple 
pole  in  \x  —  Xm\  remains.  This  quantity  is  integrable  over  the  two-dimensional  surface  of  the  tile, 
and  it  is  evaluated  analytically.  The  resulting  nonsingular  expressions  are  suitable  for  numerical 
evaluation. 


An  intermediate  result  crucial  to  obtaining  the  Laurent  series  is  given  by 


n-u 
lim  - 

X-^Xm  U 


A  ^  A  -  -  - 

-2“  ■ 


(Cl) 


where  u  =  x—Xm  and  is  the  curvature  tensor  P=  (1  —  nn)-Vn  evaluated  at  x  =  Xm-  The  third 
order  tensor  x  is  related  to  the  second  derivative  of  the  normal.  Terms  containing  Xzm  eventually 
vanish  when  angular  integrations  are  performed.  For  a  sphere  of  radius  R, 


P 


R 


(C2) 


Equation  (Cl)  was  derived  by  using  a  coordinate  system  with  the  2-axis  normal  to  the  surface  at 
Xm-  An  expansion  in  the  surf2K:e  height  was  performed,  and  finally,  the  results  were  converted  to 
a  coordinate-independent  form. 


With  Wm  as  defined  in  Appendix  B,  and  assuming  for  the  moment  that  no  air-sea  interface 
is  present,  the  diagonal  matrix  elements  are 


Br) 

Cr, 

Dr, 


(1  -  — )  Tr(p)  Wrr,  -  n(x„)  •  Vfci„t(x^) 

loTT  \  Pint  /  oTT  Pint 

f  (^eit  “  fcmt(Xm)) 

A 

47r 


(j^ext  ^inti^rn)^  (j^ext  ^fnt(^m))  Sm 

+  i  -^Tr  {yr^TkintiSm))  +  ^Tr(p)  n[xm)  ■  Vkint{xm) 
(1  -  — )  TV(p)  Wm  +  ^^niXm)  ■  VkmtiSm)  Sm. 

lOTT  \  Pext )  oTT  Peit 


(C3) 

(C4) 


(C5) 

(C6) 


The  operator  Vx  is  the  derivative  tangential  to  the  surface.  Recall  that  the  ansatz  used  to  generate 
the  interior  Green’s  function  implicitly  made  use  of  the  WKB  approximation.  In  the  version  of 
the  approximation  used,  the  second  derivative  of  kmt  was  neglected,  and  it  is  therefore  consistent 
to  do  so  in  the  equations  above. 
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Equation  (28)  is  obtained  by  combining  the  above  results  with  additional  terms  resulting 
from  the  presence  of  the  air-sea  interface.  The  latter  are  nonsingular  and  their  contributions  to 
the  surface  integral  are  found  by  simply  evaluating  the  corresponding  integrands  at  the  centers  of 
the  tiles  and  multiplying  by  the  respective  tile  areas. 
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WIREFRAME  ANNEALING 


Our  wireframes,  originally  designed  for  optimal  fit  on  a  sphere,  suffer  distortion  when  they 
are  projected  onto  nonspherical  plume  shapes.  In  this  appendix,  we  describe  the  method  used  to 
correct  that  distortion  by  adapting  the  wireframe  to  the  new  shape. 


The  goal  of  the  procedure  is  to  produce  a  nonspherical  wireframe  whose  tiles  are,  as  nearly 
as  possible,  identical  equilateral  triangles.  Generally  this  can  only  be  achieved  in  an  approximate 
way,  so  the  key  is  to  quantify  the  phrase  “as  nearly  as  possible.”  For  this,  we  first  number  the 
tile  edges  (the  line  segments  connecting  the  adjacent  vertices)  with  an  index  j  =  1,2,  •  •  • ,  L  and 
denote  the  length  of  the  jth  edge  by  ej.  The  edge  lengths  then  have  a  mean  v'alue  of 

(Dl) 

^  j=i 

and  an  absolute  deviation,  relative  to  that  mean,  of 


E  = 


7  >  0  • 


j=l 


(D2) 


We  note  that  E  reaches  its  theoretical  minimum  value  of  zero  for  the  ideal  wireframe  of  identical 
equilateral  triangles.  If  the  vertices  are  allowed  to  move  about  on  the  plume  surface,  the  configu¬ 
ration  that  minimizes  E  will  approximate  that  ideal  as  nearly  as  possible  given  the  shape  of  the 
plume  and  the  topology  of  the  vertex  connections.  Thus  we  adopt  £"  as  an  objective  function  to 
be  minimized. 


With  this  objective  function  in  hand,  we  need  an  efficient  method  for  sampling  the  vertex 
configurations  and  converging  to  the  minimum.  The  method  we  chose  is  a  novel  synthesis  of 
two  well-known  techniques  —  the  downhill  simplex  method  and  simulated  annealing  —  that  was 
recently  implemented  as  the  AMEBSA  algorithm*.  It  provides  a  natural  way  to  use  annealing 
methods  in  continuous  control  spaces.  For  this  approach,  a  wireframe  configuration  is  specified  by 
the  angular  coordinates  of  all  the  N  vertices.  Each  configuration  may  be  thought  of  as  a  point 

w  =  d2,4>2,  •  •  ■  ■,9x,<j)N)  (D3) 

in  a  control  space  of  2N  dimensions.  To  initialize  the  procedure,  a  simplex  of  such  configurations 
(a  set  of  2N  -)-  1  points  in  the  2Ar-dimensional  space)  is  constructed.  The  bare  downhill  simplex 
algorithm  would  operate  by  applying  a  set  of  “moves”  (reflections,  contractions,  expansions,  or 
multiple  contractions)  to  the  vertices  of  this  simplex,  proceeding  in  the  “downhill”  direction  defined 
by  E{yf).  AMEBSA  enhances  that  procedure  by  using  simulated  annealing  to  inject  an  element 
of  randomness  into  these  moves.  The  result  is  an  algorithm  that,  as  the  annealing  temperature  T 

♦W.  H.  Press  and  S.  A,  icukolsky.  '‘Simulated  Annealing  Optimization  over  Continuous  Spaces."  Comput.  Ph\s.. 
Jul/Aug  1991. 
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is  gradually  lowered,  takes  the  system  point  efficiently  downhill  along  any  narrow  valleys  in  the 
control  space  and  at  the  same  time  avoids  getting  stuck  in  local  minima.  The  procedure  is  aptly 
described  by  AMEBSA’s  authors  as  follows*: 

At  a  finite  value  of  T,  the  simplex  expands  to  a  scale  that  approximates  the  size 
of  the  region  that  can  be  reached  at  this  temperature,  and  then  executes  a  stochastic, 
tumbling  Brownian  motion  within  that  region,  sampling  new.  approximately  random 
points  as  it  does  so.  The  efficiency  with  which  a  region  is  explored  is  independent  of  its 
narrowness  (for  an  ellipsoidal  valley,  the  ratio  of  its  principal  axes)  and  orientation.  If 
the  temperature  is  reduced  sufficiently  slowly,  it  becomes  highly  likely  that  [the]  simplex 
will  shrink  into  that  region  containing  the  lowest  relative  minimum  encountered. 

The  effectiveness  of  the  method  for  wireframe  optimization  can  be  seen  in  the  following  tw'o 
figures  that  were  generated  during  the  evolution  of  the  annealed  ellipsoidal  wireframe  of  Fig.  7  from 
its  unannealed  predecessor  in  Fig.  6.  Figure  D1  shows  the  progress  of  the  procedure  as  a  function 
of  the  number  of  annealing  iterations.  While  the  annealing  temperature  is  reduced  monotonically 
(according  to  an  empirically  determined  cooling  schedule),  the  objective  function  E(w)  decreases 
in  a  regular  way,  “freezing”  into  its  minimum  value  after  about  100  iterations.  Figure  D2  displays 
the  areas  of  the  tiles  in  the  final  annealed  wireframe.  The  “northern  hemisphere”  tiles  are  on  the 
left  in  this  figure  and  the  “southern  hemisphere”  tiles  are  on  the  right.  The  left  and  right  halves  are 
nearly  identical  only  because  the  plume  shape  on  which  the  annealing  was  done  is  symmeiric  across 
its  equator.  The  dotted  line  marks  the  average  tile  area  and  the  dashed  line  indicates  a  theoretical 
refe^'ence  value  —  the  plume  area  divided  by  the  number  of  tiles.  The  average  is  satisfyingly  close 
to  the  reference  value  and  there  are  only  a  few  statistical  outliers.  These  outlier  tiles,  the  ones 
with  areas  above  0.01  (tile  numbers  near  1  and  100)  turn  out  to  be  the  eight  tiles  connected  to 
the  north  pole  and  the  corresponding  eight  at  the  south  pole.  The  annealing  implementation  used 
in  this  report  did  not  adjust  these  polar  tiles  effectively.  Although  we  did  not  think  the  effect 
significant  enough  to  justi^  recomputing  the  results  in  this  report,  we  have  devised  an  improved 
scheme  for  future  use. 


Fig.  D1  Relative  variance  and  temperature  as  functions  of  the  number  of  annealing 

iterations 


♦ibid. 
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COMPUTATIONAL  PRACTICALITY 


This  appendix  investigates  the  practical  issues  involved  in  the  computer  solution  of  the 
discrete  continuity  condition  as  proposed  in  Section  5  of  this  report. 


As  noted  there,  the  continuity  condition  is  a  2iV-dimensional  linear  system  of  the  form 


2Nx2N2Nxl  2Nxl 


with  a  known  constant  matrix  M.  First  b  is  calculated  and  then  the  system  is  solved  for  x.  For 
spefed  and  efficiency,  it  is  desirable  to  do  this  entirely  “in  core.”  In  that  case,  the  need  for  (2iV)^ 
complex  RAM  storage  locations  for  M  is  the  critical  factor.  This  storage  requirement  grows  with 
the  square  of  the  number  of  surface  tiles.  Since  the  linear  dimension  of  each  tile  must  remain 
smaller  than  an  interior  wavelength,  it  is  clear  that  the  number  of  them  needed  to  form  a  surface 
wireframe  will  be  proportional  to  the  square  of  the  acoustic  frequency.  Thus,  the  storage  needed 
for  M  will  increase  in  proportion  to  /'*. 


The  wireframe  approach  will  always  be  practical  at  low  enough  frequencies,  but  it  will  rapidly 
exhaust  the  RAM  available  in  any  given  computing  environment  as  higher  frequencies  are  con¬ 
sidered.  For  a  quantitative  estimate  of  the  maximum  practical  frequency,  we  consider  a  spherical 
plume  with  a  diameter  D.  If  Cjnt  is  a  typical  interior  sound  speed,  then  the  corresponding  interior 
wavelength  at  frequency  /  is  Xint  =  ^nt/f-  The  edges  of  the  tiles  can  be  no  larger  than  a  frac¬ 
tion  /3Xint  of  this  length,  where  0  ~  0.25  for  instance,  so  that  the  coherent  details  of  the  internal 
field  can  be  adequately  sampled.  The  area  of  a  t3q)ical  tile  can  be  expressed  as  a{0\int)^  where 
a  <  1.  A  square  tile,  for  example,  would  have  a  =  1.0,  while  an  equilateral  triangle  would  have 
a  =  sin(60°)/2  =  0.43  .  Thus 


(E2) 


For  the  present  numerical  estimate  we  use  triangular  tiles  whose  sides  are  one  quarter  of  an  interior 
wavelength  on  average;  i.e.,  a  =  0.5,  0  =  0.25.  Tables  El  through  E3  list  the  resulting  values  of 
2N  for  £>  values  ranging  from  1  to  3  meters  and  /  in  the  band  [50, 300]  Hz  with  Ci„t  =  100,  200, 
and  300  m/s.  If  we  suppose  that  the  available  RAM  is  8  MB,  then  at  4  bytes  per  single  precision 
real  number  (as  on  a  VAX  machine),  M,  a  single-precision  complex  matrix,  can  be  held  “in  core” 
provided  that 


2N  < 


/8  X 10® 

\  2x4  J 


=  1000. 


Inspection  of  these  tables  confirms  that  most  realistic  plume  diameters  and  frequencies  can  be 
handled. 
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Table  El  -  2N  for  typical  /  and  D  values  with  cj  =  100  m/s 


■1 

50  Hz 

100  Hz 

150  Hz 

200  Hz 

250  Hz 

300  Hz 

1  m 

50 

201 

452 

804 

1256 

1809 

2  m 

201 

804 

1809 

3216 

5026 

7238 

3  m 

452 

1809 

4071 

7238 

11309 

16286 

Table  E2  - 

2N  for  typicad  /  and  D  values  with  cj  =  200  m/s 

■■ 

50  Hz 

100  Hz 

150  Hz 

200  Hz 

250  Hz 

300  Hz 

1  m 

12 

50 

113 

201 

314 

452 

2  m 

50 

201 

452 

804 

1256 

1809 

3  m 

113 

452 

1017 

1809 

2827 

4071 

Table  E3  - 

2N  for  typical  /  and  D  values  with  q  =  300  m/s 

■■ 

50  Hz 

100  Hz 

150  Hz 

200  Hz 

250  Hz 

300  Hz 

1  m 

5 

22 

50 

89 

139 

201 

2  m 

22 

89 

201 

357 

558 

804 

3  m 

50 

201 

452 

804 

1256 

1809 
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PARTIAL  WAVE  COMPUTATION 


This  appendix  deals  with  acoustic  scattering  from  a  spherical  fluid  body  whose  interior  sound 
speed  is  a  function  of  the  radial  coordinate  only.  The  analytic  approach  to  the  problem  is  briefly 
reviewed  first  and  then  a  numerical  implementation  is  developed.  The  close  analogy  to  quantum- 
mechanical  scattering  from  a  spherically  symmetric  potential  is  exploited  by  following,  wherever 
possible,  the  analytic  treatment  of  the  problem  in  Mertzbacher’s  popular  text  [FI]  and  by  adopting 
the  notation  used  there. 

The  scattering  body  is  embedded  in  a  imiform  fluid  medium  and  both  of  them  have  the  same 
density.  The  acoustic  field  for  this  situation  obeys 

+T{r,k))u{T,k)  =  0  (FI) 

where  k  =  uifcext  is  the  exterior  wavenumber,  n{r)  =  Cext/c{r)  is  the  refractive  index,  and  T(r,  k)  = 
‘n?{r)k‘^  is  the  square  wavenumber  at  r  —  the  analog  of  the  kinetic  energy  in  the  Schrodinger 
equation. 

FI.  PARTIAL- WAVE  ANALYSIS 

In  spherical  coordinates  r  =  (r,  $,  <^),  the  partial  wave  series 

oo  +( 

Uir,  fc)  =  E  E  UimYP^ie,  ct>)ut{r,  k)/r  (F2) 

<=0  m=—t 

applies  with  each  radial  eigenfunction  Uf  (r,  k)  satisfying 

ug{r,  k)  +  Kt{r,  k)u({r,  k)  =  0  (F3) 

where 

Ke{r,k)=nr,k)-e{e  +  l)/r^  (F4) 

is  the  “kinetic  energy”  with  the  centrifugal  term  included.  When  f  >  0,  the  differential  equation 
is  singular  at  the  origin.  The  physical  solution  is  the  one  for  which  ue{r,  k)/r  remains  bounded. 

The  elastic  scattering  amplitude  has  the  partial  wave  series  F(i9,  fc)  =  where 

j?  is  the  scattering  angle  and 

Ft{'d,  k)  =  {2e  -I-  l)F<(cosi?)  X  at{k)/k  .  (F5) 

The  first  factor  is  a  known  real  quantity  and  the  second  involves  the  complex  expression 

at{k)  =  sin  6e{k)  (F6) 
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in  which  6t{k)  is  the  £th  phase  shift.  For  a  body  of  radius  a,  classical  arguments  [FI,  pages  234-235] 

dcj 

indicate  that  the  series  can  be  terminated  at  roughly  ^  =  6  =  ka. 


The  problem  of  calculating  F,  then,  reduces  to  finding  a 

finite  number  of  these  ae  terms. 

This  is  customarily  done  with  the  aid  of  the  real-valued  auxiliary  quantities  ^eik),  S({k),  Ae{k), 

0e{k)  which  satisfy 

exp{i2^e{k))  = 

hf\b) 

hf\b) 

(F7a) 

st(k)  = 

1 

(F7b) 

Mk)  = 

-“‘'•‘fg 

(F7c) 

Pe{k)  = 

f  a  dut{r,k)\ 
\ue{r,k)  dr  J 

- 1  (F7d) 

r—a 

where  h!j^\  are  spherical  Hankel  functions  and  '  indicates  differentiation  with  respect  to  the 
argument.  /3({k)  is  the  log-derivative  of  the  radial  factor  ue{r,  k)/r  at  r  =  a  and  is  the  ^3^  — »  0 
“hard  sphere”  limit  of  6^.  In  terms  of  these  quantities,  Eq.  (F6)  becomes 


Qf  =  e 


i2ii 


( 


St 

0e  -  (A^  -I-  ist) 


(F8) 


Another  quantity  of  some  interest  is  the  5-matrix  eigenvalue 

Se{k)  =  exp{i25^(fc)} 


=  exp{z2^^(A:)}  ^  1  -I-  6-^  — 


.(2)' 


(b)  ') 


^(6) 


m)  -  b 


(b) 


which  is  the  ratio  of  outgoing  and  incoming  amplitudes  in  the  ^th  radial  factor. 


(F9a) 

(F9b) 


It  is  convenient  to  formulate  the  numerical  computations  in  terms  of  ordinary  cylinder  func¬ 
tions  by  exploiting  the  following  standard  expressions  [F2| 

qwix)  =  (Fio) 

QM  =  -Q.+i{x)  + -Qu{x)  (Fll) 

X 

in  which  v  is  any  real  index,  Q  is  one  of  the  cylinder  functions  J,  F,  and  q  is  its  spherical 

counterpart  among  Specific  forms  used  in  the  following  are 

h^P{x) 
h^p{x) 

hT'jb)  _  (hf^'{b)Y 

h^^\b)  V  hf{b)  ) 


(x) 


^e+i/2\-^) 


bH. 


(1) 

<+1/2 


ib) 


(F12) 


(F13) 
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The  computation  of  proceeds  as  follows; 

•  Find  exp{i2^t{k)}  by  Eq.  (F7a), 

•  Compute  ^i{k)  =  arg  (exp{i2^£(fc)})  /2, 

•  Compute  st{k)  by  Eq.  (F7b), 

•  Evaluate  0e{k)  by  solving  Eq.  (F3), 

•  Find  Ae{k)  from  Eq.  (F7c),  and 

•  Find  a({k)  from  Eq.  (F8). 


To  aid  in  the  interpretation  of  results,  one  may  also  compute: 

•  St{k)  by  Eq.  (F9b) 

•  6^{k)  =  arg(exp{i2«£(fc)})/2. 


U{k),  st{k)  and  Ae{k)  require  nothing  more  than  some  well-known  special  functions.  The  crux  of 
the  matter  is  to  compute  f3e{k)  by  integrating  the  radial  equation  from  r  =  0  to  a. 

F2.  SOLUTION  OF  THE  RADIAL  EQUATION 

With  k  and  £  dependence  suppressed  for  the  moment,  the  definitions  v  =  u'  and 

=  ( :[:! 

=  (  J( 

put  the  radial  equation  into  the  form 

U'(r)  —  A(r)U(r 

The  evolution  of  U(r)  from  some  point  r  =  p  to  r  =  a  can  be  expressed  in  terms  of  P(a,  p),  the 
propagator  matrix  [F3] 

U(a)  =  P(o,p)U(p).  (F16) 

P(a,  p)  performs  the  transformation  from  initial  values  given  at  p  to  final  values  at  a  and  usually 
needs  to  be  found  by  numerical  integration.  For  this  we  use  an  adaptation  of  the  trusty  fourth 
order  Runge-Kutta  routine  in  Numerical  Recipes  [F4j.  The  implementation  is  fairly  straightforward 
with  special  attention  needed  for  only  two  details,  namely  the  singularity  of  the  equation  at  the 
origin  and  the  numerical  instability  that  can  occur  in  intervals  where  K{r)  is  negative.  These  are 
discussed  below. 


1 

r)  0 


1  =  0 


(F14) 


(F15) 
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F2.1  Singularity 

For  ^  =  0,  the  origin  is  a  regular  point  of  Eq.  (F3),  and  Eq.  (F16)  can  be  invoked  with  p  =  0 
using  the  initied  value 

'0\ 


U(0)  = 


1 


(F17) 


For  ^  >  0,  this  cannot  be  done  because  the  origin  is  a  singular  point.  Instead,  we  first  find  an 
analytic  approximation  for  the  regular  solution  U(p)  at  a  point  p  near  the  origin,  and  then  apply 
Eq.  (F16).  The  Taylor  expansion 


c(r)  =  c(0)  +  7r  H - , 

where  7  is  the  radial  derivative  at  the  origin,  implies  a  similar  series 

27 


nV)=n^(0)  + 


(F18) 


(F19) 


for  the  square  of  the  refractive  index.  Thus,  2|7|r/c(0)  estimates  the  fractional  error  incurred  in 
approximating  by  its  value  at  the  origin.  To  put  it  another  way,  for  a  set  tolerance  6,  the 
estimate  n^(r)  «  n^(0)  is  valid  out  to  a  radius 


r  = 


c(0)6 
2|7|  ■ 


The  approximate  radial  equation  in  the  interval  [0,f],  i.e.,  Eq.  (F3)  with 

Kt{r,k)  =  T{Q,k)-i{t+l)lT\ 
has  a  regular  solution  that  is  easily  found  to  be  proportional  to 


■) 


(F20) 


(F21) 


(F22) 


To  a  fixed  tolerance  77,  then,  u  oc  and  u/v  =  p/[t+  1)  out  to  a  radius 


4^  +  6 


yT{Q,k)  • 

Thus,  finally,  we  find  that  the  proper  initial  value  for  Eq.  (F16)  is 


(F23) 


U(p)  = 


P 

€+  1 


(F24) 


at  the  starting  point  p  =  min(f,f). 


F2.2  Stability 

Because  the  eigenvalues  of  A(r)  are  p(r)  =  ±y/-K{r),  the  character  of  the  solution  U(r) 
depends  on  the  sign  of  K{r).  When  it  is  positive,  the  fundamental  solutions  that  make  up  U(r)  are 
both  sinusoidal.  Since  neither  of  them  grows  to  dominate  the  other,  the  numerical  ode-solver  has 
no  trouble.  For  negative  K{r),  however,  the  fundamental  solutions  are  exponential,  with  one  of 
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them  growing  and  the  other  decaying.  Whatever  the  computer  wordlength,  if  K  remains  negative 
over  a  large  enough  interval,  these  exponentials  will  eventually  differ  so  greatly  that  the  ode-solver 
will  be  unstable  due  to  roundoff  error.  In  the  present  case,  the  most  convenient  remedy  for  this  is 
to  generate  the  solution  in  a  series  of  sufficiently  short  subintervals  and  “re-orthogonalize”  [F5]. 

An  estimate  of  what  is  “sufficiently  short”  can  be  based  on  the  case  where  K  is  constant 
over  an  interval  [rj,  r/]  of  length  w  =  |r/  —  r,|.  Here  the  fundamental  solutions  starting  from  ri  are 
exp(-t-|r  -  ri\n)  and  exp(— [r  —  ri|p).  For  initial  values  to  be  faithfully  propagated  to  r/  without 
significant  roundoff  problems,  the  wordlength  must  be  sufficient  for  the  computer  to  distinguish 
exp(— 2u;/i)  from  0  -I-  exp{—2wn).  The  requirement  is  exp(— =  Be  where  e  is  the  machine 
roundoff  and  B  is  a  large  number.  Thus, 

WmaxitJ-)  = -^og{Be)/2n  (F25) 

is  the  largest  interval  that  the  ode-solver  can  be  expected  to  handle  correctly.  For  e  10“®  (typical 
of  VAX  machines)  and  B  =  10^,  for  instance,  Wmaxipi)  =  3.5//X. 

In  the  general  nonuniform  case,  we  denote  the  maximum  eigenvalue  over  the  interval  by  p 
and  refine  the  interval  into  N  equal  subintervals  (ri,ri],  [ri,r2],  •  •  • ,  [rAf_i,r/],  each  much  shorter 
than  WmaxiP)-  With  that  done,  we  may  compute  the  solution  at  r/  without  significant  roundoff 
problems  by 


U(r/)  =  {P(r/, rN-i){P{rN-i,rN-2)  ■  ■  •  {P(r2, n){P(n,  J’i)}}  •  •  'DUlri)  (F26) 

in  which  the  ode-solver  is  invoked  to  compute  each  of  the  P’s.  Since  P(r,rj)  is  the  solution  to  the 
initial  value  problem  [F3] 


P'(r,rj)-A(r)P(r,rj)  =  0  (F27) 

P{rj,rj)  =  1,  (F28) 

we  can  generate  P(rj+i,rj)  by  simply  applying  the  ode-solver  twice,  with  each  application  using 
one  column  of  the  unit  matrix  as  initial  values.  The  order  of  operations  mandated  by  the  curly 
brackets  is  crucial.  Tire  or. idling 

U(r/)  =  {P(r/,  rAr-i){P(r^-i,  riv-2)  •  •  •  {P(r2,  ri){P(ri,  ri)U(ri)}}  •  ■  •}},  (F29) 

while  analytically  identical  to  Eq.  (F26),  is  numerically  equivalent  to  U(rf)  =  P(7'/,ri)U(ri)  and 
thus  does  nothing  at  all  to  alleviate  roundoff  troubles. 

F3.  SAMPLE  RESULTS 

For  the  following  results,  the  scattering  body  is  one  meter  in  radius  and  the  exterior  speed 
is  1500  m/s.  Figure  10  displays  |F(7?,  A:)!  together  with  /:)|  for  £  =  0,  ••■,3  as  functions 

of  frequency.  The  scattering  angle  is  45°  and  the  interior  sound  speed  is  a  uniform  500  m/s. 
Resonances  are  to  be  expected  near  the  frequencies  /„  =  n  x  125  Hz  indicated  by  a  “back  of 
the  envelope”  calculation  using  Cext  =  oo.  These  are  confirmed  in  that  figure.  Figure  FI  shows 

the  modulus  of  the  total  scattering  amplitude  as  a  function  of  frequency  for  the  same  situation 

except  that  the  interior  sound  speed  is  a  linear  function  of  radius.  The  average  sound  speed 
(c(0)  +  c(o))/2  is  held  at  500  m/s  as  in  Fig.  10  by  shifting  c(0)  and  c(a)  by  multiples  of  10  m/s 
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FREQUENCY  (Hz) 


Fig.  Fl  -  |F|  as  a  function  of  /  at  =  45“.  c<nt(r)  is  linear  with  Cj„t(0),Cjnt(a)  taking  on 
the  following  values;  (A)  500,  500  m/s;  (B)  490,  510  m/s;  (C)  480,  520  m/s. 


in  opposite  senses.  The  figure  illustrates  a  maxim  of  scattering  theory:  the  average  sound  speed 
determines  the  lowest  resonance  and  the  finer  details  of  Cintif)  impact  the  higher  resonances. 
Figure  F2  shows  the  dependence  of  the  total  scattering  amplitude  on  angle  as  well  as  frequency 
for  the  same  case  investigated  in  Fig.  10.  In  this  figure,  /  is  varied  in  1  Hz  steps  over  the  range 
[50  Hz,  500  Hz]  and  extends  firom  10°  (near  backscatter)  to  170°  (near  forward  scatter)  in  10° 
increments.  The  resonamces  are  clearly  visible  as  ridgelike  structures  parallel  to  the  angle  axis. 
Because  of  the  symmetry  of  the  problem,  each  of  them  can  be  associated  with  a  unique  t.  As 
expected,  the  £th  resonance  exhibits  an  (^  +  l)-Iobed  angular  structure. 


Fig.  F2  -  |F|  as  a  function  of  /  and  t?  for  the  same  case  as  Pig.  10 
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Appendix  G 

RELATION  TO  EARLIER  WORK 


This  appendix  provides  further  details  about  the  place  of  our  work,  within  the  context  of 
earlier  efforts. 


Gl.  RELATED  PROBLEMS 


Among  the  scattering  problems  that  are  close  kin  to  ours,  the  simplest  are  the  exterior 
problems.  These  deal  with  impenetrable  targets.  Some  fixed  boundary  condition  is  imposed 
on  the  surface  and  the  scattered  field  is  inferred  in  the  exterior.  The  first  of  such  problems 
to  be  solved  involved  scattering  firom  symmetric  bounded  shapes.  Reference  Gl,  for  example, 
derives  the  scattering  amplitudes  for  shapes  such  as  spheres  and  spheroids  with  either  haurd  or  soft 
boundary  conditions  as  closed-form  expressions  involving  special  functions.  But  that  approach  is 
practical  in  only  a  few  cases.  A  more  comprehensive  approach  is  described  in  Ref.  G2,  which  treats 
electromagnetic  scattering  firom  perfect  conductors  of  arbitrary  shape,  addressing  both  analytic 
methods  and  computer  implementation.  An  extension  of  this  type  of  problem  is  the  treatment  of 
scattering  firom  bossed  surfaces  [G3,  G4].  In  such  problems,  boundary  conditions  are  applied  on 
an  infinite  planar  target  with  irregularities,  the  bosses,  protruding  from  it  at  regular  or  random 
locations.  Here  again,  classical  methods  based  on  special  functions  can  produce  valuable  results 
when  the  bosses  have  simple  shapes.  Any  scattering  situation  is  essentially  an  exterior  problem  if 
the  target  is  either  totcilly  absorbing  or  totally  reflecting,  but  these  two  extremes  do  not  cover  the 
important  cases  in  which  the  target  is  penetrable. 

Scattering  from  a  penetrable  object  is  necessarily  an  interior/exterior  problem.  On  the  sur¬ 
face  of  the  object,  physical  continuity  conditions,  rather  than  extinction  conditions,  are  applied 
and  these  serve  to  connect  the  interior  and  exterior  wave  solutions.  In  such  problems,  the  target’s 
internal  structure  is  important.  Penetrable  targets  with  spherical  symmetry  have  been  analyzed 
with  great  success  using  methods  that  were  perfected  in  the  quantum  mechaniczJ  study  of  scat¬ 
tering  from  radial  potentials  [G5].  The  main  effect  of  the  interior  structure  is  the  occurrence 
of  resonances  —  analogs  of  the  nearly  bound  “virtual  states”  in  atomic  scattering.  For  individ¬ 
ual  bubbles,  the  lowest-frequency  resonance,  a  monopole  feature,  has  been  analyzed  with  special 
thoroughness  [G6].  Scattering  from  spatial  arrays  of  such  resonant  monopoles  has  received  con¬ 
siderable  attention  [G7,  G8]  and  collective  effects  have  been  recognized.  These  particular  analyses 
have  been  limited  to  small,  precisely  positioned  arrays  and  have  focused  on  coherent  phenomena 
such  as  “superresonance”  [G9-G11].  However,  incoherent  collective  effects  in  large-scale  ensembles 
of  elementary  scatterers  have  also  been  analyzed.  As  will  be  described  in  more  detail  shortly,  the 
fundamental  acoustic  result  of  introducing  large  numbers  of  bubbles  into  water  is  an  incoherent 
effect  —  the  reduction  of  the  effective  sound  speed  in  the  medium.  Scattering  from  plumes  formed 
of  such  anomalously  slow  bubbly  water  clearly  is  an  interior /exterior  problem.  Traditional  meth¬ 
ods  are  adequate  for  simple  shapes  with  homogeneous  interiors,  but  more  realistic  cases  call  for  a 
theoretical  approach  that  is  both  comprehensive  and  numerically  implementable. 
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G2.  BOUNDARY  INTEGRAL  EQUATION  METHODS 

Both  direct  and  indirect  numerical  approaches  are  possible  for  the  general  scattering  problem. 
Direct  methods  operate  by  first  generating  a  3-D  grid  throughout  the  target  and  some  of  the 
surrounding  medium  and  then  applying  a  finite  difference  solution  method  on  the  grid.  This 
straightforward  approach  has  its  drawbacks.  One  of  these  is  the  large  number  of  mesh  points 
needed  for  the  3-D  grid.  Another  is  the  strong  dependence  of  solutions  on  the  dissipation/dispersion 
properties  of  the  finite  difference  method  that  is  used.  Boundary  integral  equation  (BIE)  methods 
circumvent  those  problems  by  an  indirect  approach.  They  need  far  fewer  mesh  points  since  they 
use  a  2-D  grid  (on  the  target  surface).  They  also  readily  incorporate  the  far-field  limit.  BIE 
methods  begin  by  formulating  the  scattering  problem  in  terms  of  3-D  volume  integrals.  They  then 
exploit  Green’s  integral  identities  to  convert  these  to  2-D  integrals  over  the  target  surface.  'T’he 
main  historical  drawbacks  connected  with  the  BIE  approach  have  been  (a)  nonuniqueness  of  the 
solution  at  interior  resonance  frequencies  and  (b)  troublesome  (1/r^)  singularities  in  the  surface 
integrands.  Recently,  however,  it  has  been  discovered  that  a  more  adroit  BIE  formulation  can 
eliminate  both  of  these. 

The  earliest  applications  of  BIE  methods  were  to  exterior  scattering  problems.  As  Lamb 
explained  long  ago  [G12,  article  290],  direct  application  of  an  elementary  BIE  approach  in  this  case 
fails  to  produce  a  unique  solution  at  the  eigen-firequencies  of  the  corresponding  interior  problem. 
However,  this  no  longer  constitutes  a  barrier  to  the  use  of  the  method.  Techniques  have  been 
developed  to  overcome  the  nonuniqueness  difficulty  where  it  occurs  [G13,  G14]  so  that  BIE  methods 
are  now  in  common  use  [G15-G17].  Besides,  many  important  exterior  problems  lie  outside  the 
troublesome  frequency  range  [G18]. 

It  should  be  emphasized  that  the  exterior  scattering  problem  itself  is  mathematically  well- 
posed.  When  nonuniqueness  occurs  in  BIE  solutions,  it  is  an  artifact  of  the  particular  way  the  BIE 
method  is  formulated.  The  same  is  true  of  the  interior/exterior  scattering  problem.  Here,  however, 
techniques  to  preserve  uniqueness  have  never  been  widely  used.  Indeed,  they  have  not  been 
needed  at  all  since  Kittappa  and  Kleinman  [G19]  developed  an  interior /exterior  BIE  formulation 
with  explicitly  unique  and  stable  solutions  [G20].  Actually,  Ref.  G19  demonstrates  uniqueness 
for  only  a  restricted  range  of  wavenumbers,  but  the  proof  is  extended  in  Ref.  G21  to  all  realistic 
wavenumbers. 

Reference  G19  also  provides  a  means  of  “regularizing”  the  singular  kernels  in  the  surface 
integrals,  i.e.,  isolating  the  singular  parts  into  a  separate  term  which  is  manifestly  well-behaved. 
This  regularized  expression  appears  to  be  a  promising  starting  point  for  further  analysis.  We 
do  not  use  it  here  because  it  does  not  seem  to  facilitate  numerical  computations.  In  numerical 
applications,  the  weak  singularity  of  the  kernels  has  been  handled  either  by  boundary  displacement 
[G15]  or,  more  recently,  by  variational  reformulation  of  the  problem  [G18,  G22].  Since  neither  of 
these  techniques  seemed  appealing  from  a  practical  standpoint  for  targets  with  arbitrarily  specified 
shapes  and  interiors,  we  identified  and  dealt  with  the  singularities  using  our  own  method. 

G3.  TREATMENT  OF  BUBBLY  MEDIA 

The  rigorous  theory  of  wave  propagation  in  media  containing  large  random  distributions 
of  scatterers  (bubbly  liquids,  for  instance)  has  been  formulated  in  the  full  nonlinear  by  Caflisch 
et  al.  [G23].  The  result  essentially  validates  the  earlier  heuristic,  physically  motivated  results  of 
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van  Wijngaarten  [G24].  Prosperetti  et  al.  [G25)  have  developed  an  approach  that  is  especially 
suited  to  describing  the  internal  dynamics  of  the  bubbles  themselves.  The  linearized  theory  was 
presented  by  Fouldy  [G26]  and  somewhat  later  by  Lax  [G27].  Over  the  past  decade  and  a  half, 
Prosperetti  and  his  co-workers  have  further  developed  the  mathematical  language  needed  for  sound 
propagation  through  a  bubbly  medium  [G28-G31].  This  work  has  shown  that  a  bubbly  mixture 
can  be  modeled  as  a  fluid  characterized  by  a  complex  effective  index  of  refraction  [G30].  For 
frequencies  wp'l  below  the  resonance  frequency  of  the  individual  bubbles,  the  index  of  refraction 
is  real  and  given  by  the  well-known  Minneart  formula  [G32,  G33]. 

G4.  SCATTERING  FROM  NEAR-SURFACE  BUBBLE  STRUCTURES 

More  recently  this  group  has  applied  these  results  to  studies  of  the  related  questions  of 
ambient  noise  due  to  bubble  structures  in  the  ocean  and  acoustic  scattering  from  these  structures 
[G32,  G34,  G35].  As  a  first  calculation,  Prosperetti  et  al.  [G32]  considered  a  hemisphere  under 
a  flat  surface.  They  solved  the  problem  by  the  method  of  partial  waves  and  used  the  results  to 
estimate  backscattering  strengths  in  the  ocean,  obtaining  values  compatible  with  the  Chapman- 
Hajris  curves.  The  effects  of  radial  layering  in  such  clouds  were  also  considered  and  w’ere  found  to 
be  of  minor  importance. 

Sarkar  and  Prosperetti  [G35]  examined  the  sensitivity  of  the  scattering  cross-section  to  the 
shape  of  the  cloud  (again,  for  clouds  attached  to  the  surface).  They  first  used  a  T-matrix  method. 
This  approach  can,  in  principle,  provide  scattering  solutions  for  hemispheroids  or  cones.  However, 
it  turns  out  that  the  T-matrix  technique  is  not  very  well  behaved  numerically.  The  situation 
deteriorates  as  the  frequency  is  decreased,  and  the  problems  also  get  worse  the  greater  the  difference 
between  the  cloud’s  radius  and  depth.  To  handle  this  situation,  they  obtained  a  perturbative 
solution  valid  for  small  cloud  depths,  where  the  bubble  structure  is  approximated  by  a  boss  rather 
than  by  a  penetrable  cloud. 
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